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SUMMARY 


Convergence difficulties have been encountered in the deflection 
analysis of certain structures with most of the bending material 
in the skin using the stiffness method The difficulties are 
attributed to the use of the elementary beam spar matrix. Two 
new spar matrices are developed and compared with the matrix 
based on elementary beam concepts. The error, as a function 
f the nodal breakdown and the structural parameters, is deter- 
mined for a beam analyzed using the spar matrix for the web and 
stringer matrices for the flanges. Convergence to the final stiff 
mess coeflicients is shown to be much more rapid using either of 
he new matrices for structures with most of the bending mate 
jal in the skin 

An unswept box beam is analyzed to demonstrate convergence 
f the stiffness solutions and the superiority of the new moment 
shear spar matrix in an actual problem characteristic of a wing 
he stiffness method results are shown to converge to analytical 
esults obtained with elementary bending theory including trans- 
verse shear deformation and shear lag corrections 


SYMBOLS 


12E/,/GAL?, shear deformation parameter re- 
ferred to a spar section 
12E/,/GAI*, shear deformation parameter re- 
ferred to beam span 
area of spar cross section effective in resisting 
shear, in.? 
beam width, in 
Young’s modulus, psi 
V/5GN 2Eb?, shear lag parameter, in.~! 
shear modulus, psi 
spar depth, in 
bending moment of inertia, in. ! 
beam span, in 
spar section span, in 
number of spar sections 
1 — 5/,/61)~', shear lag parameter 
= generalized coordinate 
Q, R, S, T = stiffness coefficients 


Received March 3, 1958 

t The material given in this paper is an extension of that pre- 
Pnted in reference 1. 

* Research Engineer, Structural Dynamics Staff, Seattle Divi- 
on. 


spar thickness, in 
displacements in the x, y, and ¢s directions, reé 
spectively, in 
strain energy, in.-lb 
strain transformation matrix 
forces in the x and ¢ directions, respectively, Ib 
matrix relating stress and strain components 
shear strain 
displacement vector 
error 
normal strain 
Poisson's ratio 
normal stress, psi 
shear stress, psi 
coordinate transformation matrix 
Subscripts and Superscripts 
cross section 
elementary beam spar matrix 
beam stations 
number of beam sections 
plate 
spar 
shear-moment spar matrix 
shear-only spar matrix 
beam tip station 
directions of the x and y axes, respectively 


INTRODUCTION 


iy THE STIFFNESS METHOD, the structure is regarded 
as an assemblage of parts. Each part is assigned 
a stiffness matrix relating forces and deformations at 
its terminals, called nodes. The stiffness matrix for 
the complete connected structure is written from the 
component stiffness matrices. A given column of the 
matrix consists of a listing of the forces at each of the 
nodes for unit deformation of the given node. Where 
two or more members have a common node, forces are 
simply added. The inverse of the stiffness matrix is 
the matrix of influence -coefficients. 

For a given nodal network, the accuracy of a struc- 
tural analysis using the stiffness method depends upon 
the choice of the component stiffness matrices. Each 
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Fic. 1. Assumed deformation modes of spar section for shear- 


moment spar matrix. 


matrix is derived on the basis of a set of admissible 
modes of deformation, and, in the solution, only a 
linear combination of these modes is possible. Simple 
strain states may be used for each of the structural com- 
ponents providing an adequate representation of the 
most complicated strain variations as the structure is 
divided into finer and finer segments. 

Judgment may be used in choosing stress, strain, or 
deformation patterns which are believed to conform 
with the principal actions of a particular structural 
For example, in analyzing spar webs, 
and normal 


component. 
shear stress may be assumed constant 
stress linearly varying throughout the depth of the 
spar; no horizontal cuts need be made in the spar web. 
By introducing these assumptions an appreciable re- 
duction in computational labor is achieved at the pos- 
sible risk of some loss of accuracy in the analysis. 

In the derivation of the stiffness matrix for the tri- 
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angular plate in reference |, simple strain states hay 
Uniform shear and constant in-plane nor 
The stringer stiffness 


been used. 
mal strains have been assumed. 
matrix can be shown to be a particular case of the plat, 
matrix—i.e., the case involving normal strain in oy 
direction only and no shear deformation. The spar 
matrix of reference 1, however, has been derived on th, 
basis of elementary beam concepts. 
been used to avoid a refined nodal breakdown of th 
spar web and flanges. 

Actual experience with the stiffness matrix method 
has demonstrated that, even with fine node networks 
considerable error may occur in practical wing analygjs 


1. Use of th 


using the stiffness matrices of reference 1. 
plate matrix only in problems of in-plane deformation 
has given excellent results. Analyses of sheet-stringer 
problems using the stringer and plate matrices hay 
yielded satisfactory comparisons with experiment an 
theory. 
vergence difficulties are attributable to use of the ek 
mentary beam spar matrix. 

This paper will consider the derivation of other spar 
matrices. It will evaluate the convergence character 
istics of these spar matrices and the spar matrix based 
on elementary beam concepts for the analysis oi 
stressed-skin wings. The examples in this paper are 
constant depth structures. 
taper alone will be discussed in a subsequent paper. 


The problems peculiar t 


THE SPAR MATRIX 


The stiffness matrix discussed in reference | for the 
spars (and ribs) is derived by treating the spar as an 
elementary beam and including transverse shear effects 
The spar matrix represents two independent solution: 
of the differential equation of bending which approxi- 
mately depict the manner in which the spar will act 


in a wing. This matrix is of the form 


X ; 6EI, L3(1 + a) | (L7/h?) + [Cl + a) L?/3h?] symmetric lu, 
3 (L?/h?) — [C1 + a)L?/3h?| (L2/h?) + [((1 + @)L?/3h?| u, 
Z; (L/h) (Lh) l W, 
Z; — (L/h) — (L/h) —1 1 u 
where 
a = 12E/,/GAL? 
X, Z = forces in the x and z directions, x positive to To facilitate comparisons, this matrix has been t¢ 
the right and z positive down (forces in written from reference 1, Eq. (11a), with the coordinat 
Ibs. ) system shown in Fig. 1. Uniform shear is assumed 1! 
u,w = displacements in the x and z directions, in. the web. It should be noted that this matrix pertail' 
E = Young’s modulus to the top flange nodes of the spar. Equal and opp’ 
i, = bending moment of inertia of the spar site X forces and u displacements and equal Z force 
L = spar span and w displacements are implied at the lower flang 
h = spar depth nodes. 
G = shear modulus If more degrees of freedom were permitted, a mo! 


A = web cross-sectional area effective in resist- 
ing shear 

subscripts denoting the left and right ends 
of the beam segment, respectively, as 
shown in Figs. | and 3 


exact spar representation could be provided. Sp 
stiffnesses could be generated using plate elements ft 
the spar webs and either plate or stringer stiffness 
for the flanges eliminating some of the assumptiot 
inherent in elementary beam concepts. Such an 4 
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EVALU 





ATION 


proach would involve considerably more calculations 


than is generally warranted in most analyses. 

The moment of inertia of the spars of some recently 
designed high-performance wings has comprised only a 
small part of the total bending material. Also, trans- 
verse Shear effects are generally small and need be ap- 
proximated to only a first order of accuracy. Therefore, 
it seems reasonable to limit spars to vertical cuts only 
and to account for shear by assuming it to act uni- 
jormly over the cross section. 

Recognizing that an exact spar matrix simulating 
point loadings is neither desired nor required, a matrix 


based upon the above precepts can be derived. A sim- 

ple spar matrix developed by superimposing the strain 
X GA/2L | (Lh?) + (aL? 3h?) 
X; (L? h?) — (aL? 3h? 
Zz, (L/h) 
Zz; —(Lh) 


This matrix will be referred to as the shear-moment 
spar matrix as opposed to the elementary beam spar 
matrix of Eq. (1 

The aflinity between Eqs. (1) and (2) is evident if 
their application is considered as a function of network 
fineness alone. As the network becomes finer, | a be- 
comes much less than 1, and second-order terms may 
be neglected in expanding {1 + (1/a)]~!' in a binomial 
Substituting this approximation in Eq. (1) 


series. 
vields Eq. (2) identically. 

In an actual analysis, it is desirable to represent long 
spans satisfactorily for a coarse network. In this case, 
the distinction between the elementary beam spar 
matrix and the shear-moment matrix is that bending 
is the primary action in the former and shearing the 
primary action in the latter. 


A third spar matrix to be denoted the “‘shear-only”’ 
spar matrix can be obtained from Eq. (2) easily by 
letting the bending stiffness vanish (a = 0). Using the 
assumptions of the previous paragraph, the same shear- 
only matrix is derivable from Eq. (1). 


PLATE-SPAR COMPATIBILITY 


In application of the method, stiffness solutions 
using the elementary beam spar matrix have been 
found to be as much as 100 per cent too flexible when 
compared to calculations based on conventional ‘‘elas- 
tic axis’ methods. One structure having 90 per cent of 
its bending stiffness in the skin and 10 per cent in the 
spars gave this poor correlation despite a relatively fine 
breakdown (75 nodes). This increased flexibility has 
been accounted for in the poor compatibility between 
the spar and plate in the interval between the end nodes. 

If the elementary beam spar matrix is considered as 
primarily a bending portrayal, the spar-plate system 
mathematically acts as shown in Fig. 2. Since the 
plates are given no out-of-plane stiffness and since they 
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Fic. 2. Incompatibility of plate and elementary beam spar 
states associated with uniform shear and a constant 
moment and considering rigid translation and rotation 
is derived in the Appendix. Fig. 1 shows the assumed 
deformation modes. The resulting spar stiffness ma- 


trix is expressed as 


symmetric u 

(L2 h?) + (aL? 3h?) u 
(2) 

(Lh | Ww 

—(L h —]| l Ww 


are only developed by being stretched between nodes, 
they may be considered plane surfaces between nodes, 
each plate cut acting as a hinged line. 

Essentially, the vertical bending stiffness of the com- 
posite structure is obtained from the plate acting as 
spar flange material. Using the elementary beam 
spar, if the moment of inertia of the spar is much less 
than the effective moment of inertia of that section of 
the wing, the spar can assume large curvatures with 
relatively little extension of the chord distance between 
nodes. In this case, a large number of spar-plate inter- 
connections are required to fully represent the com- 
patibility between plate and spar. 

To analyze the spar-plate compatibility problem, 
consider the structure shown in Fig. 3. The problem is 
to determine the effect of increasing the connections 
between the spar and stringers running along its top 
and bottom flanges. This approximates the spar-plate 
interconnection problem when a sufficient number of 
nodes have been used so that strain is sensibly con- 
stant over an interval between nodes. 

To obtain a solution, stiffnesses will be summed for 
and B. 
To make answers comparable, the rows and columns 
associated with nodes between A and B will be reduced 


an increasing number of sections between A 


out. These reductions imply there will be no imposed 
loading except at A and B. 

The stiffness matrix for a stringer is of the well- 
known form 


bo ; 1 —1 u; 
i f-ee[ i ae 


where .!.. = cross-sectional area of the stringer. 


If the stringer stiffness matrix of Eq. (3) is added to 
the elementary beam spar matrix, Eq. (1), the shear- 


(2), or the shear-only spar matrix 


\=—/» 


moment matrix, Eq. 
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for a section 7—j, the combined stiffness matrices can 
each be represented as 


bY Q symmetric U; 
X; R Q u, 
= (4) 
Zi S S T Ww; 
A —-S -S -T T W; 


/ 
where Q, R, S, and 7 are stiffness coefficients dependent 
upon the choice of spar stiffness matrix. 

Combining two such matrices to describe the stiff- 
ness relations for the nodes 7 and k and the common 
node j, and then reducing out the rows and columns 
associated with node /, yields again a four by four matrix 


Xi QO — (1/2)[(R?/Q) + (S?/T)| 

X, —(1/2)[(R2/O0) + (S?/T)]  -O — (12) [(R? 
Zi (1/2)S[{1 — (R/Q)] (1/2)S[{1 — 
Zy —(1/2)S[1 — (R/Q)] —(1/2)S[1 — 


where the subscript k denotes beam station k. 
By comparing Eqs. (4) and (5), recursion relation- 


ships are evident: 


Qon = Qn — (1/2)[(Rn?/Qn) + (Sn?/T,) | 

Roy = —(1/2)[(Ru®/Qn) + (Su/Tr)] Cg 
Son = (1/2)S,[1 — (Rn/Qn) |] | 

Ton = (1/2)[T, — (Sn?/Qn)] 


where the subscript » denotes the number of beam 
segments. 

Substituting for the stiffness coefficients in Eqs. (6), 
the stiffnesses can be found for any number of seg- 
ments. For example, the diagonal stiffness coefficients 
in a Z row using the elementary beam spar matrix are 


TT.” = (6EI,/2L*) }(1 + a) + 
[3/(1 + I,/1,) |f - 
Ti = (GEI,/4L%) }(1 + a) + 
(15/0 + 7,/1) If ™ 
Te = (GEI,/SL*) {(1 + a) + 4) 
(63/(1 + 7,/T.) If 
T, = (6EI,/nL*) \(1 + a) + 


[(n? — 1)/(1 + I,/T,) |} 


where /, = moment of inertia of the stringers about 
the beam neutral axis neglecting the moment of the 
stringers about their own centroids, and eb as a super- 
script identifies the elementary beam spar matrix. 

Writing Eq. (7) in terms of the parameters of the 
beam as a whole, L = //n and a = na, and including 
expressions obtained in a similar way for the shear- 
moment and shear-only matrices, gives 


T,@ = (6EI,/l®) {a + (U1 + 1,/n7l,) + 
(1 + J,/7,)|f? 
T,!" = (6EI,/l*) {a + [(1 — 1/n?) + 
a+1,/1)yj-{ © 
T,” = (6EI,/I) {a + ((1 — 1/n?) + 


(I,/1.)- | 


where the superscript sm identifies the shear-moment 
spar and so the shear-only spar matrix. 
Examination of Eqs. (8) reveals that when n becomes 


SPAC 
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infinite the elementary beam spar matrix and the shear- 
moment matrix yield identical results, as expected 
For values of /, /, typical of a wing it can be seen that 
the shear-only matrix gives results within a few per cent 
of those of the other spar matrices for 7 infinite. These 
conclusions can be verified for the other stiffness coefi- 
cients. 

Finally, if the error using m sections in the analysis is 
defined as A, = (T. T,,)/T., the error relations 
to Eqs. (8) are 


corresponding 


A,” = (1/n*) ;a[(./Ip) + 1] + 
(1,/Ip) + ( n*)i | 
A,” = —(1/n*) jal) + U,/7 |) + 9 
1 — (1 n*)} | 
A,” = —(1/n’) [a4U,/1,) + 1 — (1/n?)] 


Assuming n finite, a comparison of Eqs. (9) exposes 
the relative merits of the three spar matrices. Sine: 
the denominator of the shear-moment error is always 
greater than that of the shear-only, the shear-momen! 
spar matrix has superior convergence characteristics 
The absolute value of the error of the eb matrix and thi | 
sm matrix can also be compared. It can be seen that 
if /,/[, is less than one, the denominator of the sm erro: 
is always greater than that of the eb matrix and, agai! 
the shear-moment matrix exhibits superior convergent: 
characteristics. 

In Fig. 4, the error is plotted using the elementar' | 
beam and shear-moment spar matrices for variow 
In plotting, transverse shee! 

This graph makes apparel! 


values of » and J,/J,. 
rigidity is infinite (a = 0). 
the large errors that can result using the elementar 
beam spar matrix for the analysis of wings with mos 
of the bending material in the skin. 

The sign of the error functions is significant. Tl 
positive error function of the elementary beam sp 
matrix indicates that the stiffness coefficients will | 
crease with an increasing number of nodes. Experiemt’| 
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with the plate stiffness matrix in problems of in-plane 
deformation has shown that the plates have a negative 
error function. Using the plate matrix and the ele- 


Te) 
oO 


mentary beam spar matrix together can result in the 
stiffness coefficients varying as shown in Fig. 5 as the 
number of nodes is increased. If solutions were per- 


@ 
re) 


formed in the region of A and B sections, the method 
would yield an erroneous impression of convergence. 


N 
oO 





This manner of convergence was reported in the ex- 
ample of reference 1. The negative error function of 





o 
oO 


the shear-moment matrix is desirable for use with the 
plate matrix because it leads to monotonic convergence. 





uo 
oO 





ANALYSIS OF A Box BEAM 








A box beam of aspect ratio 7 is selected to demon- 
strate the convergence of the stiffness method and the 





superiority of the new shear-moment spar stiffness 
matrix for practical wing structures. This box beam 


O 
Oo 
| 


Mh 

oO 

| 
poe a 


is similar to the example used in reference | except the 





: : : > 4 . - —+ 4 
ratio of /, /, has been increased from 2 to 64. The . 


small amount of spar material was chosen to accentuate VN | ; 4 
. ; be s 
the spar-plate interaction problem, and it is somewhat 


A.,=ERROR IN 22 STIFFNESS 
ee ae ae 
(eo) 











re) 
y 





Tp | 


ratio wings being designed today. The beam is shown Is{_ | | ‘ae 


in Fig. 6. 12345678910 
The properties assumed for the beam are: 
N-NUMBER OF BEAM SECTIONS 


representative of many of the very thin, low-aspect 











Material: aluminum, & = 10.5 X 108 psi, G = nt 
3.94 & 108 psi Fic. 4. Error in stiffness coefficients for elementary beam and 
aes ps shear-moment spar matrices 

Spars: /, = 64.314 in.', 1 = 2.544 in.? 


Plate: /, = 4116.151 in.* 

Total Beam: J = 4180.465 in.! 

The box beam selected for study accentuates the 
second-order bending effects of shear lag and trans- 


verse shear. The large unsupported chord is associated EB 
with large shear lag effects while the small amount of T 
web material was chosen to emphasize transverse shear nN 


deformation. The beam is assumed to deflect without , | 
plate buckling. | | 

An elementary bending theory solution for the de- | | 
flection of the box beam consists of Reissner’s shear lag Al 








correction” and a modification for transverse shear de- 


formation. The equation for the tip deflection of a NUMBER OF SPANWISE SECTIONS 


uniform beam due to the load Z, as the tip is 


Fic. 5. Possible convergence of vertical stiffness coefficient with 


j iE — 1 sinh 4 / elementary beam spar matrix. 


ZEIT pel f cosh fl} GAS 
Samer shear lag) + Y 


where the parameters .V and / are defined as 





Pg X 1Z 








N = [1 — (51, 61)|-', f = (1b) V5NG/2E 4=400 
and 6 is the beam width in inches. h -9544" FF SPAR 
| Substituting numerical values into Eq. (10) and let- j 
ting the load be two Ibs. (i.e., one Ib. at the tip of each ey 
spar), the deflection becomes w, = 0.00127 in. Of the 
total deflection, elementary bending accounts for 76.3 Fic. 6. Geometry of box beam 


per cent, shear lag for 17.4 per cent, and transverse 
Shear for 6.3 per cent. 
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Fic. 7. Nodal breakdowns for analysis of box beam 
Unfortunately, the closed solution given by Eq. (10) 
is subject to several restrictions. It applies to a beam 
unswept untapered plan 
The problem 


Further, 


of rectangular cross section, 
form, and uniform spanwise properties. 
selected conforms with these restrictions. 
Reissner restricts the cover skin to act as an ortho- 
tropic plate. These orthotropic assumptions, described 
by Hildebrand in reference 3, are that the elastic modu- 
lus /, in the transverse direction is infinite and Pois- 
son's ratio associated with strains in the longitudinal 
direction is zero. The stress-strain relationships re- 
duce to 


é¢- = 0,/E,, ¢& =0, Yxuy Fay lG (11) 


where ¢ = normal strain, o = normal stress, y = shear 
strain, 7 = shear stress, and the subscripts refer to the 
coordinate system. 

To make the stiffness method solution equivalent, it 
was therefore necessary to alter the two-dimensional 
stress distribution developed in the plates and constrain 
the action of the general plate stiffness matrix to ortho- 
tropic behavior. This is readily accomplished. 

Three nodal breakdowns were chosen to demonstrate 
stiffness method; 15 nodes, 28 
These breakdowns are illustrated 
the stiffness 


convergence of the 
nodes, and 45 nodes. 
in Fig. 7. For each nodal breakdown, 
matrix and tip deflections were calculated first with the 
elementary beam spar matrix (including transverse 
shear) and second with the shear-moment spar matrix. 

The tip deflection resulting from these calculations 
is plotted in Fig. 8 against the number of spanwise cuts. 
A curve is drawn through the points for illustrative 
purposes. Convergence of the elementary beam spar 
matrix results is shown to be exceedingly slow for a box 
beam with this low a percentage of spar bending mate- 
rial. For 15 nodes, the tip deflection is 72 per cent 
Even at 45 nodes, a fairly fine breakdown, the 
Since the 


high. 
box beam solution is 18 per cent too flexible. 
plate stiffnesses are overestimated, the increased flexi- 
bility due to lack of compatibility between the spar and 
the plate between nodes is actually greater than the de- 
flection error. 

The tip deflection determined with the shear-moment 
spar matrix shows rapid convergence and excellent 
agreement with the analytical solution. At 15 nodes, 


the stiffness method solution is 6 per cent too stiff, while 
at 45 nodes, the solution is only 1.5 per cent stiff. 


This 
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SCIENCES SEPTE 
discrepancy is attributed primarily to lack of conyerg 
ence of the plates. 

In each of the breakdowns chosen, the ratio of chord 
wise to spanwise sections is the same namely, tyy 
Therefore, the error can be expressed as a function oj 
the number of spanwise sections only. Extrapolating 
the beam solutions on the assumption that the error js 
inversely proportional to the square of the number oj | 
spanwise sections used, in accordance with Eq. (9), th 
elementary beam solution yields 0.00127 in the limit 
and the shear-moment 0.00126. Thus, this problen 
also demonstrates convergence to essentially the exact 
solution in the limit for either spar matrix. | 

It is of interest to note the increased sensitivity t 


error of the matrix in inversion as the number of nodes 
is increased. This sensitivity may be measured by the 
ratio of the maximum to minimum eigenvalues of th, 
matrix. Values of the ratio for two by two matrices 
of 15, 28, and 45 node solutions are 1.533, 1.735, 1.846. | 
and 1.997, 2.013, 2.017 for the elementary beam and | 
shear-moment solutions, respectively. Both solutions 
will have the same eigenvalue ratios in the limit. Th 
use of this or an equivalent ratio as a single measure 
for convergence evaluation is suggested. 
stiffness method was 


In addition, the also used 


to determine the difference between the deflection 
resulting from the assumption of orthotropic plates and 
that from the more general plate action which includes 
transverse normal strains. The deflections allowing 
transverse normal strains were from | to 3 
for each breakdown and each spar matrix. 
tively, it might be anticipated that the more general 
action would yield the more flexible solution. The 
deflections of the beam with orthotropic plates are in 


Poisson's ratio effect and are 


per cent less 
Qualita 


creased by the lack of 
decreased by the rigidity in the transverse direction 
The net effect for this problem is a very slight increase 
in overall deflections using the orthotropic plates. It 
is observed that when the stiffness method is used, the | 
amount of effort required to obtain the more general 
solution is no greater than that under the orthotropic 
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EVALUATION OF 


CONCLUSIONS 


The elementary beam spar matrix and the shear- 
moment matrix yield the same results in the limit as- 
suming similar shear distribution in the web. 

The shear-moment matrix has superior convergence 
characteristics as well as giving the least absolute error 
jor a given breakdown for wing analyses where the cover 
skins furnish the majority of bending stiffness. 

Extrapolation of stiffness solutions proved feasible for 
this simple problem and merits study for practical 


wings. 
APPENDIX 


A spar matrix is desired with nodes in the planes of 
the skins only for a structure untapered in depth. 
Deformations will be constrained to equal and opposite 
u and equal w displacements for top and bottom nodes 
at each end of the spar. These assumptions are con- 
sistent with structure possessing a horizontal plane of 


V = [£t/2(1 — v?) 





l 


SS leer | |v 


0 


[Et 2(1 — v?) 





where v is Poisson's ratio, |’ is total energy, and ¢ is 
spar thickness. | is obtained by an integration over 
the spar volume. Here, the first integration has been 
performed assuming stress constant through the thick- 
ness of the spar. 

Now, using Eqs. (A-1) and (A-2) in Eqs. (A-3), the 


energy appears in the matrix form: 


V = [Et/2(1 — v2)| X 
SS (6¢-"W'V¥We- 8,) dx dz 





(A-4) 


where 


V= (Gt2L)f Sf 6,’ F [8 + v)s?2/Lh?] + (Lh?) 
—[S(1 + v)s?/Lh?] + (L/h?) 
L/h 


—(1/h) 


[S(1 


where G = E/2(1 + »). 

Integrating (A-5) and forming the partial derivatives 
of the strain energy with respect to the node displace- 
ments yields the node forces in accordance with Castig- 
liano’s theorem. The resulting matrix defining the 
relationship between node forces and displacements is 
denoted the shear-moment stiffness matrix and appears 
as Eq. (2). 

The forces are obtained for only the upper spar 
nodes. The complete system of node forces is shown 
in Fig. 3. Vertical forces are divided equally between 
the upper and lower spar nodes and the X forces in the 
upper skin plane have corresponding equal but opposite 
forces in the lower skin plane. 


SS het + Wee + €2 + (1 2) 


SPAR MATRICES 
Four modes of deformation satisfving the 
Mode 


is rigid rotation. 


symmetry. 
constraint requirements are shown in Fig. | 
Mode (b) 


and (d) involve, respectively, 


(a) is rigid translation. 
The elastic modes (c 
deformation under uniform shear and deformation due 
to a pure couple. 

The assumed admissible strain states and the corre 


sponding deformations are 


€ 0 0 O —2s L qi 
é, 00 0 22 L2 || % i \-1 
q 
7 0 0 I/L 0 gq 
u 0 (12h 0 0 ] 1 
Mu, 0 (12h O ALI | @ | 
W; 1 —(12)L 0 O q \-) 
w, 1 (2h 1 1 Ys | 
or, 6, = gq; and gq; = ¢~'6 


where g; is a generalized coordinate. 
Eqs. (A-1) and (A-2) assume the condition of plane 
stress. The strain energy can be expressed 


[1 — vly,.7{ — 





y 0 le, \ (4-3) 
l 0 . dx an 
0 (1/2)(1 — »)] Eve 
l v 0 
Yy =e | 0 
0 O (1/2)(1 — pv) 
The inverse of ¢ is discovered to be 
L/h 0 1 O 
1 | 2/h 0 0 0 
a ~ | —(L/h) -—(L/h) -1 1 
— (L/h) Lih 0 O 


Substituting the required matrices in Eq. (A-4) and 
performing the matrix multiplication gives 


symmetric 6; dx dz 
+ v)s?/Lh?] + (L/h?*) - 
(A-3) 
Lh LL 
—(1/h) —(1/L) W/L 


REFERENCES 


' Turner, M. J., Clough, R. W., Martin, H. C., and Topp, L. J., 
Stiffness and Deflection Analysis of Complex Structures, Journal 
of the Aeronautical Sciences, Vol. 23, No. 9, p. 805, September, 
1956. 

? Reissner, E., Analysis of Shear Lag in Box Beams by the Prin- 
ciple of Minimum Potential Energy, Quarterly of Applied Mathe 
matics, Vol. 4, No. 3, p. 268, October, 1946 

> Hildebrand, Francis B., Zhe Exact 
Problems in Flat Panels and Box Beams 
Transverse Direction, NACA TN 894, 1943 


Solution Shear-Lag 
Assumed Rigid in the 


of 





A Three-Dimens 
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in Subsonic and Supersonic Flows' 
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SUMMARY 


A small-perturbation theory for the description of the steady 
state flow fields of axial compressor blade rows is developed by 
superposition of cylindrical wave functions. The relation be 
tween the full three-dimensional solution and the usual strip- 
theory (cascade) approximation for blade rows having a large 
number of blades is determined, and a similarity law for three 
dimensional correction effects is provided. Numerical data 
necessary for the application of the general solution to configu 
rations of arbitrary hub ratio and large blade number are pre- 
sented, and this information is used in computing numerical ex 


amples in subsonic and supersonic cases 


SYMBOLS 


speed of sound in undisturbed fluid 


a = 

A = B/3(1 — h) = similarity parameter 

B = number of blades in blade row 

c(r = local blade chord 

Cat? = axial projection of local blade chord 

by = coefficients for determining A, 

( = pressure coefficient referred to dynamic pres- 
sure of axial flow, Eq. (55) 

C, = pressure coefficient according to modified 
strip theory 

for. x = source distribution function 

o(x = profile shape function 

h = ry/rr = hub ratio 

V2 = modified Bessel function of first kind, order v 

Ss = Bessel function of first kind, order 1 

| = radial eigenvalues 

k = summation index 

L = blade spacing, Eq. (30) 

m = summation index 

Ma: = axial Mach Number 

M(r = relative Mach Number at 7 

My = relative Mach Number at tips 

n = mB 

N, = Neumann function (J), in reference 5), 
order n 

p = static pressure 

?.. = static pressure of undisturbed fluid 
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g(r = strength of radial source lines 

q = (u, ?7,, %) = perturbation velocity vector 

OmB—i A = Fourier expansion coeflicients 

Onn-x" = numbers related to Fourier cvoeflicients, F 
(29 

ru = radius of hub 

ry = radius of shroud = radius of tips 

? = radius of sonic cylinder 

R = cylinder functions, order 7 

t = f, = time 

l = axial advanced speed 

Ai = discontinuity in velocity component normal 
to blade surfaces 

(x1, "1, 0 = cylindrical coordinates in space-fixed syster 

(es, £0) = evlindrical coordinates in rotor-fixed systen 

¥ = r70 

a, B = variables in asymptotics of evlinder functions 

B = (1 — M,,?)'/" 

5, = phase shift angles 

r, = axial a parameter, Eqs. (7 

bu = (M, | dee 

wl? = WEG yl 

By? = (M,? — 1)!/? 

&. ¢ = variables in asymptoties of evlinder functions 

p = density 

p = density of undisturbed fluid 

o = wl l 

v7 = wl] l 

o = wr,/U = B/M, 

r(r = local blade thickness-to-chord ratio 

T = thickness-to-chord ratio at hub 

eX, 7%, 4,4) = perturbation velocity potential in spac 
fixed system 

g(x, 7, 0 = perturbation velocity potential in rotor-fixe 
system 

?1 = stagger angle at tips 

w = rotor angular speed 

( )ep = cascade theory 

a = three-dimensional correction to — modifi¢ 
strip theory 

t ) = quantity referred to coordinates aligne 


with relative velocity at tips 


INTRODUCTION 


t bow IS THE FIRST of two papers describing a three 
dimensional theory for predicting the flow fielé 
induced by an axial compressor blade row of arbitrary 
blade shape and hub ratio, operating over a wide rangt 
of axial and relative Mach Numbers. In this analysis 
special emphasis has been put upon three-dimensional 
effects not predictable by ordinary cascade strip theory 
and, in particular, upon the transonic blade row, [ot 
which strip theory fails entirely. 

In the present paper, the general method will be set 
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RELATIVE WIND 


| Blade row geometry and notation, showing evlindrical 


FG 
Helical arrow indicates incoming 


coordinates fixed in) rotor 
flow viewed from rotor 


forth and its applications to subsonic and supersonic 
cases described. Discussion of the application of the 
theory to transonic cases will be reserved for the second 
article. 

It has been traditional for many years to study axial 
compressor blade rows operating in either the subsonic 
or the supersonic regime within a quasi-three-dimen- 
sional, or strip-theory, approximation. According to 
this approximation, the detailed effects of the blade 
properties can be predicted by applying theoretical or 
experimental cascade results at each radial station; in 
other words, one assumes two-dimensional flow in 
each compressor annulus. The universal acceptance of 
this approach is a natural outcome of its outstanding 
success in providing a rational method for designing 
and predicting the performance of axial compressor 
blade rows, particularly in the subsonic regime. 

Nevertheless, no clear procedure for estimating the 
relative importance of three-dimensional effects not 
included within the above strip theory has been avail- 
able. In particular, the variation with Mach Number 
of these effects has recently become a matter of no little 
importance, especially since compressors are now being 
designed to operate at high subsonic or low supersonic 
Mach Numbers. It is believed that the study reported 
here will provide useful information along these lines; 
indeed, a fundamental aim of the present paper will be 
to make clear the relation between the complete, three- 
dimensional theory and the strip theory that is so 


widely used. 


THE GENERAL THEORY 


Assumptions and the Differential Equation 


In the following analysis we shall consider a single 
blade row operating in an infinitely long cylindrical 
annulus (see Fig. 1). We shall assume that the com 
pressor geometry is such that only small velocities are 
induced by the blade row, and, consequently, a linear 
ized analysis of the flow field is appropriate. Accord- 
ingly, none of the higher order effects of shock waves, 
such as entropy gradients and vorticity, will appear 
in the theory. We shall further assume that the com- 
pressor moves in a perfect gas of low viscosity. The 
flow Reynolds Number will be taken to be sufficiently 
large that the vorticity introduced by shear at solid 
boundaries is confined to boundary layers of small 
thickness. Consequently, everywhere outside of these 
thin boundary layers the perturbation velocity, q 
(4, v,, %), will be written as the gradient of a scalar 


potential function —1.¢., 
g = Voe(X1, "1, %, fh) 1) 


and the usual boundary-layer assumptions will be made 

in obtaining the pressures on the boundaries. 
Application of these assumptions to the equations 

of nonviscous fluid motion yields the familiar acoustic 


wave equation in space-fixed coordinates 
V,76 = (1/a*)d,,1 ? 


where a is the speed of sound in the undisturbed fluid. 
In a cylindrical coordinate system advancing (at a 
speed l’) and rotating (at a speed w) with the com- 
pressor blades (see Fig. 1), the flow is steady, and 


Eq. (2) becomes 


(l — Moar*)err + (1 r?) [1 — (wr? a?) le — 


(2wl’ a*)gyp + Or + (1 Ne¢ Q (3) 


with A/,, = U ‘a, and 


Q(NX1, Vs, 6, i) = glX, Fa 0) 


Y= MT Uh, g = 6; —- wl; }) 


Particular Solutions 


This differential equation is not unfamiliar in fluid 
mechanics; it has previously been studied by Rott,’ 
Busemann,* and Davidson.‘ Particular solutions of 


Eq. (3) can readily be found in the form 


en. = Ay exp (ind + (w U) } (in, (1 — My,?)] + Ant X) Ral[Kas(r rr) (5) 


Here the radial functions, R,[A,,.(r rr) |, are normalized linear combinations of Bessel and Neumann functions of the 


form 


r COS 6.7, [Ay(r rr)| — sin 6.4. NalAn(r rr) Z,[Ku(r rr) | 


ay i = 


My? = (1 2) 4 Z,27(K yx) [1 — (2? Ky?) 


Where h is the hub ratio ry ‘rr. 


Mans 


— WZ, 7(hK yx) (1 — (n? 7K x?) }f 


The phase angles 6,, and eigenvalues A,, are determined from the two equations 
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R, (WK yx.) = R,(Kix) = O (G6) 


which guarantee the vanishing of the radial velocity at 
the shroud and at the hub. 


The quantities A,, appearing in Eq. (5) can be 
written 

Ane = (n/Bor) V (Kyx7/n?) — (er?/o;7) (7) 

where o7 =orr U, 6, =B Miz, B? = 1 — Maz”. a, is 


i.e., that cylinder on 
which the resultant Mach Number is 
unity, written in nondimensional form. 


the radius of the sonic cylinder 
undisturbed 
That 1s, on 


the sonic cylinder, 


o, =wr,/U = B/Maz (S) 


The Axial Dependence 

The choice of sign associated with the functions 
A», in Eq. (5) is made differently upstream and down- 
stream of a specified “‘matching plane’’ so as to satisfy 
the conditions of bounded perturbations infinitely far 
upstream and downstream of the blade row. When a 
given \,, is real, the choice is obvious; when a given 
d,. is imaginary, and 1/,, < 1, the Doppler effect must 
be considered. In terms of the steady-flow system of 
Eq. (3), the Doppler eifect requires that disturbances 
propagating upstream in the helical wave forms given 
by Eq. (5) must have a finer helical pitch than the cor- 
responding disturbances propagating downstream. !: 

In the completely supersonic case, 1/,, > 1, on the 
other hand, one simply specifies that no disturbances 
occur upstream of the blade row; both families of solu- 


tions are included downstream. ! 


Physical Meaning of the Solutions 


Once the choice of sign is made, the physical inter- 
pretation of the particular solutions (5) becomes clear. 
Depending upon whether the axial velocity component 
or the potential is required to be continuous at the 
matching plane (say the plane » = 0), each particular 
solution is either the flow due to a sheet of doublets 
(creating a discontinuity in potential) or the flow due 
to a sheet of sources (creating a discontinuity in axial 
velocity). Use of elementary flows arising from source 
sheets provides a meais for treating thickness—-i.e., 
nonlifting—problems; use of flows due to doublet 
sheets provides a method for handling the lifting case 
(see, for example, reference +). In the two cases the 
techniques applied are essentially the same. 

In the present paper, and in reference 1, only the 
thickness problem is treated in detail; the lifting case 
is left to future research for several reasons. First, 
exposition of the principles involved is simpler for the 


nonlifting problem; second, much of the numerical and 
theoretical information obtained during the analysis of 
the nonlifting problem can be carried over directly to 
the lifting case; and finally, the thickness problem is of 
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interest in its own right, particularly in the transonj, 
and supersonic regimes when wave drag associated wit} 
thickness appears. 

In the present paper, we require continuity of th, 
potential for each particular solution, and what follows 
will then be limited to application to nonlifting blag 
rows. Discussion of the extension to lifting problems 
is included in reference |. 


Superposition of the Particular Solutions 


Because of the linearity of Eq. (3), the general soly 
tion to a given boundary value problem can be ob 
tained by means of appropriate superposition of th, 
particular solutions. The method used here has bee 
to superpose the solutions in such a way as to creat 
B radial 
These radial source filaments can then be distributed 


source-like singularities of strength g(r 
over B helical sheets approximating the compressor 
blades, with a source-strength distribution function 
I(r, x) say, related to any specified blade thickness. 
One begins by writing the expression for the local dis 
continuity in mass flow across the matching plane for 


each particular solution—i.e., each source sheet: 


Am. = Aj ol U + whan & poB?Au, ) 
and if the matching plane is chosen at x = 0, 
Aten, = —An(20/U)rA ne R, [Kx (r/rr) (10 


These relations provide complete sets of orthogonal 
functions over both the 7 and 6 coordinates of the plane 
x = 0; hence, any desired distribution of mass-injection 
(or removal) at this plane can be provided by Fourier 


analysis.* In this case we specify 


Am = 0 0 ~ 2ka B 
Am = « 0 = 2ka/B: 
ae a ae Bye) (0) 
e(2kr B)+e 
r Amd6 = p.q(r) 
J (2k B)—e 


so that the mass is introduced in a manner correspond 
ing to B equally spaced radial source lines of strength 
qr’). 

At any radius the mass-injection distribution defined 
by Eq. (11) can be written as a Fourier series of the 
: for each n the function g(’ 


me 


orthogonal functions e 
can be expanded in terms of the orthogonal functions 
R, [Kur rr) 
verge everywhere except at the singularities. 
way the coefficients A,, in the formula for the potential 


The resulting double series will con 
In this 


can be determined for an arbitrary distribution 0 
source strength over the B source filaments. 

Once the expression for the flow due to these singu- 
larities is obtained, the singularities themselves can bé 
distributed in the classical manner over B helical sheets 
approximating the blades so as to represent the give 
thickness distribution of the actual compressor blading 
In the process, if c(7) is the local blade chord, the souret 

* For incompressible flow, this procedure was first carried oul 
by Ackeret* for two-dimensional cascades and by Meyer® fet 


three-dimensional configurations 
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strength g(7) is replaced by 


U?) + 1 dx: 


fir, x) V (wr? 
0 <x < c(r)/V (w'r?/U*) + I 


and with this normalization it is not difficult to show! 
that 


f(r, x) = Av, (7, X) (12) 


where Av, is the local discontinuity in velocity normal 


to each helical sheet approximating the blades. (It 


1) mas One (&) 
a me 


m=1 k=1 \mB— 


w | B clr) Ne 
in — &|¢ dé - zs | 
f trB- Jo 


where o=aor LU. 
distribution and are obtained through the formulas 


*! 


(1 or) V,(7, X 
oe h 


be 


0, B-A\N) = 


Av, (r, x) (Wor2+ 104) = 


COMPRE 


oo 
} 


) V o* + l R, [Ane (Lr rr)|d(r rn) | 


SSOR BLADE ROWS 
is merely a matter of convenience to write the speci 
fication of the profile shape, and thereby the functions 
f(r, x) and Av,(r, x), in terms of the axial position in 
stead of the position measured along the helical blades 


themselves. ) 


The General Solution 


When this procedure is carried out, the complete 
three-dimensional solution to the linearized flow of an 
axial compressor blade row is obtained in the form, 


where /,, < 1 anda = mB, 


; r j. we imB w 
Riel Kmp_i exp yimB , = — s+ a 
rr l os” l 


Ko, ") exp ) — Koj = x- i d& (13) 
'r ' l f 
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The coefficients Q,,4-,() are the Fourier coefficients of the expansion of the radial source-strength 


In the case of supersonic axial flow (.1/,, > |) we find instead 


g(r, 0, x) = 0, 


Be 
¢(r, 6, x = | » Bs 
trp" J0 m= AmB h 


exp 


where the Q,,n_,(€) are again given by Eq. (14) and 


o< ss Ic(r) Vo? + Hes} 


VN = C(r) (16) 


Yo + Thies, 
x > [c(r) Vor + Tooe| 


PROPERTIES OF THE SOLUTIONS 


The solutions given above have two particularly 
striking features that must be given close attention 
before trying to apply them to specific cases. The first 
of these is that the order of the radial functions is pro- 
portional to the number of blades in the blade row. 
Thus, for a large number of blades, we find ourselves 
seeking solutions to Eqs. (6) that involve Bessel func- 
tions of very large order. In that case, advantage can 
be taken of the asymptotic formulas developed by 
Watson,® and this procedure leads to numerical tech- 
niques that can be greatly simplified through the use 
of a correlating parameter involving the hub ratio and 
the number of blades. In fact, as we shall see subse- 
quently, the eigenvalues A,,, can be written, for large n, 


Ky n? = 1+ (2C,/n?”), (n = mB) (17) 


and in the limit of very large 1 the coefficients C;, become 
functions of Am” ” alone, where 


A = B11 — h) (18) 


i, B-} (g )R., B [ ee Br" rr) 


(14) 
Onp (AX) RipiAme (Yr) | 
(x < QO) 
Soo 7 i. 
\? sinh X,,8—, UV (x — E) x (15) 
¢ mB w 
Sime (6 n =) + mee oba, (wo 
C;~ 


As we shall see, extraction of the subsonic and super 
sonic cascade theories from the three-dimensional re 
sults is accomplished by taking the limit B ~ © with 
A (rather than /) fixed, and relaxation of this limit 
are pro- 


three-dimensional corrections that 


(2/3) 


leads to 
portional to B 
The second striking feature, namely, the possibility 
that one or more of the functions A,, may vanish |see 
Eq. (7)], is a kind of acoustic resonance and requires 
special consideration. In order to give a complete 
picture of the flow for general blade configurations, a 
method for handling such resonant cases is needed. 


However, resonance in this form is strictly a transonic 


phenomenon, since if A,, = 0, this implies 
or/os = Kun (19) 
and it can be shown’ that A,, » > 1 forO0 <u < oa 


Thus, resonance occurs only when the radius of the 
sonic cylinder is less than that of the shroud (and is 
really important only when it is actually in the flow 
consequently to 
Since application of the present 


annulus). Resonance is peculiar 
transonic operation. 
theory to transonic blade rows is the subject of the 
second of these papers, the method for treating tran- 


sonic resonance will be discussed there. 
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Fic. 2. Behavior of Bessel and Neumann functions of large 
order. Note that the first inflection point of J, is very near 


> = n and that the first zeros of both functions oceur for s > n 
Also shown is a typical linear combination of Bessel and Neumann 
functions with The phase shift is 
positive when R, leads J, 


a positive phase shift, 6,. 


SOLUTION OF THE RADIAL EIGENVALUE EQUATIONS 


The solutions (13) and (15) to the general boundary 
value problem for nonlifting blade rows are of use only 
if sufficient data are provided for determination of the 
eigenvalues A,,,, and phase angles 6, These quantities 
are determined through Eqs. (6). 

The minimum value of 7 that interests us (aside from 
n = 0, which we need not discuss here) is equal to the 
number of blades in the device under consideration. 
For axial compressors, this number is usually large, say 
on the order of 40 to 60. Because of this we can con 
veniently make use of asymptotic expressions for cylin 


der functions of large order. ! 


Properties of Large-Order Cylinder Functions 


The asymptotic formulas for large-order Bessel and 
Neumann functions given in reference 5 are diilerent 
for different ranges of the argument with respect to 
the order.* The range of argument that interests us 
most is that in which the first several zeros and ex 
trema occur, and it can be shown’ that this range is 
given for J/,(z) and V,,(s) by 

zg—n| = O(n'") (20) 
It follows that the range is the same for any linear com- 
bination of /,,(z) and .V,(z), e.g., R,(s). 

Both the Bessel and Neumann functions are mono 
tonically increasing functions for values of the argument 
between zero and the order n. An immediate conse 
quence of this is that a linear combination of these two 
functions cannot have more than one extremum in this 
range, a fact that is important in determining the region 
of possible solutions of (6). A schematic illustration of 
a Bessel and a Neumann function for large m is pro- 
vided in Fig. 2. 

* More recently, unified expressions for the asymptotic be 
havior of large-order cylinder functions have been given by Olver” 


in terms of Airy Integrals. These expressions cover the whole 


range of argument and embrace all the results listed separately 


in references | and 5. 


SCIENCES SEPTEMBER, 1958 


Dependence on the Hub Ratio, h 


The range of argument of R,, in which we are inter 
ested for any value of k is between /A,, and A Eq 
(20) implies that A, = O(n) if k is not too large (0; 
the order of itself). It will be seen that Eq. (17 
agrees with these considerations and with Eq. (20 
The fact that, at most, one extremum can occur j 
R,,(s) for s <n means that /A, 
than » but that A, must be greater than » in any case 
Consequently, the coefficients C 
of order unity. 

Now, if / is any fixed number less than unity, th 
product hk, 


to 7—1.e., 
n ) oe ai | 


if h does not involve n (or B). With the asymptoti 
formulas available in Watson's treatise’ for this rang 
(n —2>n'”), it is quite easy to show that in this cas 
the phase angles 6,,,-) are proportional asymptoticalh 
to 
” lili , 22 
Then if # is considered to be independent of B (for pur 
poses of considering the limit of large B) the phasi 
angles 6,,,-, are exponentially negligible and the func 
tions X,,,z are, for all practical purposes, pure 
functions of large order over the range considered. 
Unfortunately, pure Bessel functions of large order 
because of their property that they do not differ appre 
ciably from zero until 7 O(n'*), 3 <n, provide 
an extremely slowly converging series when used i1 
Eqs. (13) and 
computation over the whole flow annulus their ust 


would be out of the question. 


The Parameter A 


suggests a way to avoid this situation an 


Eq. (22 


vastly improve the convergence properties of the 


Fourier-Bessel series. In practical cases /i is reasonabh 


near unity. If, then, we regard (1 h) in Eq. (22 


being proportional to B for B large, we perform: 


transiormation of the argument of the radial functions 


such that the phase angles 6,,-; become of order unit) 


I. 


or at least appreciable, and the values of R,,4(A,,,,—/ 


also become significant. In other words, if we treat 
the phase angles 6,,z_; as functions of the parameter 4 
lgiven by Eq. (18)| and of B rather than functions 0 


h and B, we define a transformation of the argument 0! 


the radial functions such that the entire flow annulus 


is included within the range of argument over which the 
Ther 


the radialiunctions R,,|A,,.(7 rr) | can be used effectivel 


large order cylinder functions are appreciable. 


in representing a given distribution over the range 0 
radii with which we are concerned, and the series 1! 


Eqs. (13) and (15) will converge rapidly. 


Asymptotic Formulas 


When we treat A as a quantity of order unity, we find 


that AK,,n-, also lies within the range given by Eq 


may or may not be less 


are positive numbers 


differs from ” by an amount proportional 


Jessel 


(15), and for the purposes of practical 
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”)), In seeking solutions to Eqs. (6) we may accord- 
ingly use the asymptotic forms given by Watson? for 


this range [ypical of these expressions, we have, for 


o P n sech a 


~ [(tanh a) 3] X 


a gy [/ 1/3) (€) — iF 3(£) (23) 


3) tanh’ a 


plus errors that are bounded by terms that are specified 
in Table Here the functions /,(x) are the modified 
Bessel functions of the first kind. Watson also gives, 
for 2 n,” 

z J,(n sec 8) = [(tan B) 3] X 
m3 — 8B) -— Ff [J-as(O) + Sis | 


Sit a ti 
tan 3) 3] cos }n(tan B — B) — ct x 

J -) — Jy 73(6)]; (= (n 3) tan*® B 
[hese expressions cannot be differentiated directly 
to obtain the correct asymptotic formulas for the deriv 
itives ol so the fundamental identities relating 


derivatives of Bessel and Neumann functions to the 
functions themselves should be used in Eqs. (6). In 
each case, the error obtained by using Eqs. (23) and 
24) and their related formulas is at most of order 1 7. 
Watson provides asymptotic formulas of this type 


y 


for J,(¢) in 
n, but for .V,,(z 


this range with ¢ greater than or less than 


only for s greater than ”. In order to 


solve Eqs. (6) we also need the asyinptotic expression 


for V,,(s) when s < ”, and this has been worked out in 
reference using techniques similar to Watson's. A 
summary of the asymptotic expressions for J/,(s) and 
\V,(z) and the associated errors is provided in Table 1. 


Numerical Solutions of Eqs. ‘6 


Using the tundamental identities for the derivatives 
of the evlinder functions and expressions such as (23 
ind (24), Eqs. (6) can be put into a form in which they 


can readily be solved, with ” as a parameter, on any 


COMPRESSOR 


BLADE ROWS 549 


Eqs. (6) appear long and complicated and will not be 


reproduced here. But when x = mB is really large, 


expressions such as (23) and (24) can be simplified by 


neglecting terms of order ¢" "n ~~’ or &? "n ’, and 
then Eqs. (6) take the form, for example, when /A,,, 


mB, 


I _ (2/3) (&) — To 3(& I _ (23) (ER) — Jo36 s 
» 
FAY | 
I » ay(E + [5 2(é I _« > + Jo c 
where é, = (1/3) [2Am™” — (36,)° F (26 
and a similar expression holds for /A,,, mB! 


These equations do not involve n = mB as a parameter. 


Through them the eigenvalues A,,,z-; can be obtained 
once and for all, for sufficiently large n. These equa- 
tions can easily be solved using a desk computer and 
the tables of reference 6,f and the cor flicients C, can 


be obtained directly from ¢, through the relation 
C.(A) = (1 2) (8)° 


The more exact equations replacing Eqs. (25) and 
(26) provide small but significant corrections to the C;, 
making them not functions of Am” ” alone, as they are 
when Eqs. (25) and (26) hold, but also weakly depend 
ent on mB itself. The complete equations have been 
solved on an IBM 650 Magnetic Drum Calculator at 
the Cornell University Computing Center for three 
The results are given 
for B = 64 


1, 2, 3, along with solutions to Eqs. (25 


moderately large values of mB. 
in Fig. 3 in terms of C, plotted versus Am 
and mm and 
26) and their counterparts, which are designated by 
af ©''in the Figure. The coefficients C, determine 
the A, B 


The phase angles, obtained in the same v 


through Eq. (17 
ay, are 
plotted in Fig. 4 in terms of 4/3 tan 6 versus C,(A 
Values of the Coefficients 

Examples of values of the coefficients €) B &) tor 


a specified blade geometry were also computed by 


high-speed computer. Written in this form, however, integrating the complete relations such as (23) and (24 
Expressions (23) and (24) can be unified either using Olver’s + The corresponding numerical solutions using the unified 
results or the more recent work of Chester, Friedman, and expressions of Olver can be obtained using tables of Airy In 
Ursell I treatment given here follows that used in reference tegrals An alternative procedure to the one presented in this 
( completed before Olver’s results were brought to paper would be to develop an asymptotic series for the A, in the 
wthor ttention form of equation (1.2), page 329, reference 10 
TABLE | 
Asymptotic Formulas: nm =O 
Region The Functions ./ JN Error Bound: ] Refer 
ence 
‘ 1, (n sech a) ~ (tanh a/3) exp {| (tanh a a) +} X 3/n) exp jm (tanh a a) \0 Watson, 
a ; I é I é p 250 
, Vin sech a) ~ —(tanh a V3) exp {nla tanh a) — &} &X (O@./n) exp {nla tanh a x 
ee ; I f Th, 3(é 18 $ sinh [(7/2) tanh* a] X 
exp n/3) tanh® a 7 tanh? a 
, J,(n see B 1/3) tan 6 cos }n(tan B 8 oc} xX 24 O3/n Watson, 
= s ; 2 bol ‘ » . » 252 
ae Jats O + SJv3(O)] + C1 V3) tan 8 sin {n (tan B 8 ri > Pp. 202 
J ¢) J ¢) 
, V,.(n see 8) 1/3) tan 6 sin jm (tan 8 B ci; xX 24 0,/n Watson, 
= . 9) { ’ LD BA) ] 
} 3 J Oo+s f)| — (1/V 3) tan 8 cos |m (tan B B) —¢} X Pp. <v< 
J 4) ef 4 
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Fic. 3. Numerical information for determining the radial eigenvalues Kx [solutions of Eq. (6)]. The coefficients Cy, depend primaril) 
upon the parameter 4 = B*/s(1 — A) and secondarily on n = mB. For all n > 192 the ‘‘n = ©” curves can be used. The A,, art 
determined from the values of C;, through Eqs. (17). 
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Fic. 4. Numerical information for determination of the phase shifts, 6,4. Note that the 6,; are all positive rhe solid lines give 
values of \ 3 tan 6,,; the dashed lines give values of its reciprocal For n > 192 and C; > 8.7, use the relation tan 6,, = tan {(1/3) X 
2C r/4 
for B = 64and m = 1, 2, 3 and by integrating approxt1- the profile shape. Assuming that A 1.6 (h 0.9) 
mate expressions neglecting terms of order ¢" “n~ or oan and 
s m > 4. ; 
. . > — (r) = tolArr 1”) (h- _ ) > (28 
In this sample calculation of the coefficients, the ; en ith V Wer I) (o 
profiles were assumed to be similar and the axial pro- the numbers 
iection of the chord, c(v) V ao? + 1 = Car, to be a con- ‘ aa/j 
. +r . . . dn p—e (XxX) = Onp_w(x 1(X) To hoy? 4 (29 
stant. Then the radial and axial dependence of the Om ) Oma—e(X) lor g(x) toh h*or ! 
were obtained and are listed in Table 2. For m > 4, 








blade thickness can be separated, and 
Av, (r, x) = r(r)g(x) (27) 


where r(r) is the local thickness ratio and g(x) specifies 


TABLE 


But if 
.* given in the Table for 


these values depend only on the choice of A. 
B # G64, the values of Q,.% 
m < 3 must be modified somewhat even for A = 1.6. 


») 


Values of Qng_«* for B = 64, A = 1.6 


V h?o7r? + 1/0? + 1 


Tir) = Tf hry } 
l 2 3 4 
” 

l 0.3243 0.1089 -) OO4O 
) 4 0. 2659 0.2952 0.0341 0. OO87 
1). 2337 0.2218 Q.1111 0.0235 
0). 2131 0.2089 0.1756 0.1078 
1973 0). 1873 (0.1784 0.1441 
t 0.1853 0.1753 0.1640 0.1585 
7 0.1757 0.1658 0.1454 0 1560 
S 0.1679 0.1580 0.1384 0.1349 
7 0.16138 0.1515 0.1326 0.1294 
10) 0. 1556 0.1459 0.1275 0.1245 


5 6 az i 

0.0011 

~() O24] 0.0012 0 OO12 

0.0729 —() OO27 0.0131 0.0111 
0.1172 —() O457 0.0255 0.0168 
0.1356 0 O898 0.0229 0 0276 
0.1374 0.1247 0.06038 0 OO77 
0.1354 0.1280 0.1037 0 0293 
0.1185 0.1215 0.1181 0.0774 
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All the coefficients, of course, depend on the blade 


geometry assumed. 


APPLICATION OF THE THEORY 


The numerical information provided in Table 2 and 
Figs. 3 and 4 is sufficient to enable us to apply the 
solutions (13) or (15) to blade rows operating at spe- 
cific axial and resultant Mach Numbers in the subsonic 
or supersonic regime. Once the profile shape, and 
thereby g(x), is specified, the values of the potential, 
the velocities, and the pressure distribution can be ob- 
tained by straightforward calculation on machines such 
as the IBM 650 Magnetic Drum Calculator. This is 
actually the simplest way to proceed, from the point of 
view of obtaining numerical results. The series are 
rapidly convergent in the subsonic case and good 
accuracy is obtained very easily. In the supersonic 
case use of special techniques for series with oscillating 
partial sums improves the accuracy obtainable. 

Nevertheless, it is useful and instructive to investi- 
gate analytically the relation between the solutions (133) 
and (15) and the familiar (two-dimensional) cascade 
theory that has been so successful in predicting sub- 
sonic and supersonic compressor performance. We 
shall see that such an investigation will provide us with 
expressions for three-dimensional terms not predict- 
able by simple blade-element or strip theory. 


The Two-Dimensional Limit 


It is clear that any two-dimensional approximation 
to annular flows involves the limit of large cylinder 
radi. The relation between three-dimensional blade 
rows and two-dimensional cascades is established by de- 
fining this limit such that as r7, ~ ©, B > © simul- 


taneously, so that we have 


lim (2ar7/B) = L (30) 


rT & 


Bo « 
where L is the spacing between blades in the infinite 
cascade. In addition to this requirement, however, it 
is necessary to place a restriction on the hub ratio, and 
it is shown in reference | that, in order to obtain the 


Examples of the Two-Dimensional Limit 


Incompressible Case 


the solution (13) for .\/,, = 0. In this case, 


AmB = (1 


and we may write (where real part is implied) 


g(r, 8, x) = 


K mB—k 


tar 0 


If we put r79 = y, we obtain from Eq. (30) 


Bo = 


SPACE 


or) Kips 


Bar (Car)mar _~ ‘ im (¢) ; r é we 
= ag | Zz 2. ‘ . XG Ring (Ann k Jexp JimB (0 — ) x—€ 
2 0 m=1k=1 'r { l 'r f 


Bary ((ardmar 2 Gor(é ae 3 - bere . 
"| - 2 %, (Ku “\ exp Ku bag (35 
e rT , 


2ry i 3 


‘ 
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SCIENCES 


two-dimensional limit from the complete three-dimep. 
sional solution, one must also require 


hm hk = 1 3] 


B= 


which of course implies 7 ry — | within the floy 
annulus. 

Moreover, for reasons of convergence already dis. 
cussed, taking the limit B ~ © of expressions (1:3) and 
(15) will not be successful unless (1 — h) approaches 
°) Any slower approach will cause th 


zero as BU 
series to converge too slowly; any faster approach will 
destroy the completeness of the set of R,. For thes 
(30) should be taken 
Syvmboli 


reasons the limit defined by Eq. 
with the parameter A, rather than /, fixed. 
cally, we may describe the process as 


¢3p(B, h) ———> gon + GB, A) 32 
iq? 
B— 
1 fixed 


in preference to 
¢sp(B, h) —— Cum + @ (B, h) 33 
ry 
B— « 
h fixed 


In Eq. (32) ¢2» will correspond to cascade theory, where 
this exists, and ¢(B, A) will represent correction terms 
that are nonvanishing only upon relaxation o! the 
limit B = o. On the other hand, ¢;;,, = 
in general, and ¢)(/) does not always vanish as B > 


gop + gilh 


with / fixed. 

Taking the limit B ~ © of solutions (13) or (15) in 
the manner thus defined leads immediately to the 
classical first-order cascade theories at r = rr appro- 
priate to the Mach Number regimes involved. From 
Eq. (13) the correct cascade theory for subsonic flow, 
including the Prandtl-Géthert rule for relating the re- 
sults to equivalent incompressible flows is obtained. 
From Eq. (15) the familiar Ackeret theory can be ex- 
tracted, and, in that case, no Mach waves propagate 
forward of the cascade, since /,, > 1. The corre- 
sponding case of Ackeret theory in which the resultant 
Mach Number is supersonic (but not near unity) while 
the axial Mach Number is subsonic, so that Mach waves 
propagate forward of the cascade, is obtainable from 
Eq. (13). 


A simple example of the procedure involved in this limit-taking is provided by considering 


ie BI 


ba} Ko, 





where | 


large B 


since tl 


obtain 


since ¢C 
pressio 
velocit 
angle « 
Cc mpre 
holds f 
The 
rigorou 
missibi 
and int 
sult re 
exercisi 
inclusi¢ 
or of 
puted ¢ 
to illus 
and thi 
case. 
Subs 
dimens 
similar 
ference 
throug’ 
vert tl 
form, \ 
Mach ! 


and th 





men 


flow 


dis- 
and 
ches 
the 
will 
hes¢ 
iken 
boli- 


‘ms 
the 


4 


AXIAL COMPRESSOR BLADE ROWS 553 
B(wt/ U) = (27/L)é tan ¢7 (37 Kinp- ~ |m\B: m0 
where @r is the stagger angle at the tips and tan @7 = independently of k. Then 
Keeping only the first term in Eq. (17) for very : ‘ 
ci ss j | ; (Kmp—x/tr) |x — &| ~ (29 !m//L) {a £ 3S 
large B, we have 
Since in the limit we also have r rr @ 1, Eq. (35) be 
comes 
olrr, 8 —> o(y,x) = - DD Dd [67Onn()Rnw(K me x 
Boa J0 ” 1 / l 
exp ) (2am i) [7(y — £tan @r) — ix -€é ‘ ae B . ot Oo. (é Ro( Kor) o7 ( , ke ; 
- 1g Zz exp | —(Ku/rr) |x — &\5 d& (39 
24m lr Jo bel Ko 
If we now calculate, for example, the axial velocity component forx < 0. we have 
1 ( ens , j2am A ; ps 
v(x, 1 = pd > lor7Qmnp—n(E) Rine( Kime) | exp - [((x — &) + ty — €& tan ¢r) dé - 
i 7 0 m l1k=1 { 
= os > , 
Dd larQo(é)Ro(Ko,) | dé (40) 
re iB e ft) j l 
since the ratios Ao, rr all vanish in the limit. The sums over k can be evaluated immediately from Eq. (14), and we 
obtain 
gr(x, y) = —(1/L) Av, (rr, §) ; p eo stan 7) 4+ (1 2)¢ (dé cos gr) 
Jo ai f 
(41 
<a 
= (1/2L) } Av, (rr, —) coth } (mr L) [(x £) + i(y — Etan @r)]} (dE ‘cos 7 
/7 WW) 
since cos @r = 1 Vor? + 1. This is exactly the ex- 


pression for the axial component of the perturbation 
velocity associated with an infinite cascade of stagger 
angle @7, spacing L, and arbitrary profile shape in in- 
compressible flow, and the last expression in Eq. (41) 
holds for any x. 

The derivation given here is not, of course, fully 
rigorous. Careful attention should be given to per- 
missibility of taking the various limits inside the sums 
and integrals, and for this purpose the reader may con- 
sult reference 1. In fact, considerable care must be 
exercised before the correction terms arising from the 
inclusion of more terms in the expansion of A,,z—) 
in the compressible case) can be com- 
The present calculation merely serves 


(or of Amp 
puted correctly. 
to illustrate the relation between the three-dimensional 
and the two-dimensional theories in the incompressible 
case. 

Case—The of taking the two- 
dimensional limit in the subsonic case (\/,, # 0) is 
The only dif- 


Subson te process 
similar to the calculation given above. 
ference is that extra Mach Number terms are carried 
throughout; for example, we may use Eqs. (17) to con- 
vert the formulas for the quantities A,,z-, into the 


form, where Mp = M,, V or? + 1 is the resultant tip 
Mach Number, 
Mune = (n/Bor) V[(1 — Mr?) ‘B2] + (2C,/n* *)) 19 
(47) 
(n = mB) 
and then expand each X\,, in powers of n”~ ~~’ in the 


same way as the A, were expanded in the incompres 
sible case. When only the first term of each of these 
expansions is retained, the result is the correct subsonic 
cascade theory at r = rr, which can be reduced to the 
form (41) through the Prandtl-Géthert transformation. 


Supersonic Case—When the tip Mach Number is 
sufficiently greater than unity, both Eqs. (14) and (15 
yield the appropriate Ackeret theory in the limit de 
fined by Eqs. (30) and (31). The difference between 
the two cases, of course, is in the fact that for the one 
Mi. < 1 and for the other \/,, > 1. Inthe correspond 
ing Ackeret theory for the ./,, 
family of Mach waves propagates forward of the 


1 case [Eq. (13) 
one 
cascade; in the case given by Eq. (15) both families re- 
main behind the cascade. 

For the purposes of an example, let us consider the 


limit of Eq. (15). If we define 


= M,,’ 


ar? = Mr’? -— 1, 2’ 


we see that both ur and yu are real in this case and that 
all the \,,, are imaginary. The leading term in the ex 
pansion of these quantities is given by 


Ane ™ (—inur por) (44 
for the first several k. Substituting in Eq. (15) and 
differentiating, we have in analogy with Eq. (41), for 


oo. 
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A similar expression for the limit of the tangential velocity component can easily be found; putting the two t 
gether, the formula for the component of the disturbance velocity aligned with the resultant free-stream velocity 


vector can be constructed: | 


| *7(x) dé M,2* sin or + ur COs or ; 
u’ = ¢, cos dr + ¢,sin dr = — (| Av, (rr, §) | = — sin gr | X 
Qurl 0 cos or be 
; j2mim ' M7 tan dr + py .. i 
ye exp * (y — & tan ¢r) + - ie ==) ( 
m=— a ' L Me 
(x) dé ur cos dr — AM," sin or : 2rim 
| Av, (rr, §&) | r = + sin or z. exp ; (y — & tan ¢7) 
JV cos r Ko m = L 


( M,,* tan dr — wy ) E ) 
= (x — &) (46 
ys { 


In this case, since the arguments of the exponentials in the sums are imaginary, the integration can be performed 
directly, using the closure property of orthogonal functions. That is, if the functions y;(r) are orthonormal over 
the set /, then the closure property is written 


| fr) Skdk(r) dr = f(r’) 
JVE h 


‘ P ini ke te 
where ( ) denotes complex conjugate. In particular, if y,(r) = e* V 27, then 


(1 2x) | fim) >& 2-2) dey = f(z) (47 
x k=— 


If z = 0, + 27, ..., we have by symmetry 
| "(Qe te) te, (€20| .. =e 
f® > e™\—| dé = fleie)]; | = 0, +1, +2,... (48 
2m J t(2lx- w) k=- dé 


and in this form the closure property can be applied directly to the integrals in Eq. (46). If we identify 


Qa M,,” tan or + pw 
2) (£) = F (y — tan gr) — = = (x — &) | 
‘ " (49 
ae Qn : M,,* tan @r — ur F 
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' & we 
dzo"! 2a Mar? tan ¢r + ur 
- dé = —tan dr + = = dé 

dé Z, be F 
then ‘ : (00 

dz) 27 M,,;? tan or — ur| ,. 

di = tan or — : dé 

| dé E pe 
and the integrals in Eq. (46) become 
1 ( py : en” >¥(x) ae WE ) es 
ua’ (x, ¥) = — | Av,(rr,) 2) es | d& + } Av,(rr,) 2 ems ‘ de i” 
i 4 THT le 0 m= — dé /J0 ” dé 

, I Av,,(77, E; rr) I Av, (7, E; * ~ 
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LMT j 2 kr] ~ 
where the ranges of j and / are such that if 0 < &" < y(x) andO < &,") < y(x). &)"? and é, are defined such that 
(27/L) i (y — & tan gr) — [(M,,? tan dr + ur) wu?) (x — ~ yi = + 2j7; ¢=). 1. 2. ) = 
(0 

(29/L) \(y — &, tan or) — [(M,,? tan @r — mur) /‘u?) (x — £,2)| = +2/n; i = 0, 1, 2, f 
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in accordance with Eqs. (48) and (49). 
In terms of coordinates (x’, y’) oriented with the 


free stream, One obtains! 


£ = cos ori (x’ — ury’) + JL sin or ¥ 
urjL cos or} | (54) 
a ury’) + 
iL sin ¢r ¥ upjl cos or} | 


{0 
= cos Or) (x 


The lines (x’ + wry’) = const. are the familiar Mach 


lines of the first and second families. The extra terms 


drawn from the field point) 
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in Eqs. (54) account for the spacing between blades 
that contribute to the velocity at a given point. 
Expressions (54) used in Eq. (52) show the identity 
of the present result with the Ackeret theory, recalling 
the nonlifting nature of the blade row. Because they 
are nonlifting, the blades are what Busemann has called 
“limp bodies’’; under these conditions waves striking a 
given blade are canceled out rather than reflected, and 
the camber line of each blade is determined by the inter- 
ference from its neighbors so that there will be no load. 
Thus, Eq. (52) can be written, in accord with Ackeret, 


‘y') = —(aMyp ur) (SLOPE at upper surface of first airfoil intersected by forward Mach line of first family 


+ (aMy ur) (SLOPE at lower surface of first airfoil intersected by forward Mach line of second family 


drawn from the field point) 


and the SLOPE is given by the linear combination of 
the slope due to thickness and the slope of the induced 


camber line. 


Srrip THEORY AND THREE- DIMENSIONAL CORRECTIONS 


The above calculations correspond to assuming B = 
= ©,r rp = 1, andh = 1, without any consider- 
If, on the other 


2. ry 
ation of the approach to such a limit. 
hand, we consider that B is large but finite and that 
r ry (and /) differ from unity by an amount proportional 
to AB’ '~”’, we obtain from the full three-dimensional 
solution the theory 
known as strip theory, plus 
can be shown rigorously to be proportional to B 

We have already noted that the best way (numert- 
cally) to obtain results from the three-dimensional solu- 
tion is by straightforward calculation on a high-speed 
The best way to obtain (again 


quasi-three-dimensional 
terms that 


usual 
“correction” 


2 
2/3) 


computing machine. 
numerically) the three-dimensional corrections to strip 
theory is, consequently, to conipute the complete three- 
dimensional results and to compare these with the pre- 
dictions of strip theory. Such comparisons are pre- 
sented at the end of this paper. 


Analytic Results 


Analytic study of the relation between the strip 
theory and the full three-dimensional solution, how- 
ever, is of more than academic interest. First, it 
puts the first-order strip theory for the first time on a 
rigorous foundation and makes clear under what condi- 
tions it can and cannot be used (for example, it cannot 
be used in the transonic regime). Secondly, it tells us 
the order of magnitude of the three-dimensional cor- 
rection terms, including their approximate variation 
with Mach Number, and provides an interesting simi- 
larity law relating these correction terms between dif- 
ferent machines. Finally, it suggests an alternative 
“strip theory,’’ somewhat simpler than the usual one, 
which might in some cases give results as accurate as 
those of the ordinary strip theory with much less effort. 

The rigorous derivation of the strip-theory-plus- 
corrections expansion of Eqs. (13) and (15) is lengthy 


and need not be reproduced here. Briefly, however, 
one begins by using Eqs. (42) at any tip Mach Number 
not too near unity and expanding the first several 
Ane in powers of n°”. The leading terms of these 
expansions used in Eq. (13) or Eq. (15) with 7 rr # 1 
yield the quasi-three-dimensional strip theory, corre- 
sponding to the results given above with 7 rr; treated 
as a parameter. Keeping higher-order terms in the 
expansion of Eq. (42) and including these in Eq. (13) 
or Eq. (15) leads to the three-dimensional corrections 
previously mentioned. (Note that the condition that 
My not be too near unity implies that the necessary 
expansion of Eq. (42) cannot be done in the transonic 
regime, and therefore the strip theory has no meaning 
for transonic operation. ) 

When the first two terms of each expansion for the 
first several \,, are inserted in Eqs. (13) and (15) and 
the series over k are broken up into appropriate partial 
sums,! one may use the convergence properties of the 
series to show that, for example, 


oy 


C,(r, 9, x) = [(b — peo)/(1/2)p.U? 
—2[(u + ov,)/U] = C,,(r6, x; 7) + 


A » 


C,(r, 9, x) + of ee al (55) 
where C,, corresponds to (but is net exactly equal to) 
the usual strip-theory result and is of order t+ 

V1 —M,? C, is the three-dimensional correction 
term of order 7)B~‘* *; all neglected terms are of order 
7»B~' or smaller. The variation of C, with V7 is some- 
what complicated, and depends, for example, on the 
blade profile as well as the region of the flow field being 
For the most part, however, correction 
' when the transonic 


investigated. ! 
terms grow as Br~? = (1 — M,? 
regime is approached from the subsonic side and as 
ur? = (M7? — 1)! when approached from the super- 


sonic side. 


A Modified Strip-Theory Rule 


It has already been mentioned that C,, does not quite 
correspond to the usual strip theory, and it can be 
shown quite easily that the difference is proportional 
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to AB” *: that is, it is of the same order as e.. In 
the usual strip theory (which might be called ‘‘intui- 
tive’ strip theory) one treats y as a parameter, and the 
two-dimensional flow variables are referred to the local 
relative Mach Number and local blade profile. C,,, on 
the other hand, corresponds to a treatment that takes 
into account the blade profile at a given radius, but 
refers the two-dimensional flow field in each annulus to 
the relative Mach Number at the tips. The simplifica- 
tion obtained, particularly when the profiles are similar 
and experimental cascade results are being applied, 
could be quite significant. It is interesting that the 
three-dimensional theory leads more naturally to the 
latter rule than to the intuitive one. In the remainder 
of this paper this rule will be referred to as the ‘‘modi- 
fied”’ strip theory. 

Since the difference between the two rules is of the 
same order as the ‘correction’? terms that take into 
account additional three-dimensional effects, an inter- 
esting question arises—that is, which rule actually is 
closer to the full three-dimensional solution? This 
question can only be answered by considering the nu- 
merical results reported below. 


A Similarity Law 

Before we proceed to the consideration of numerical 
examples, however, an important theoretical result 
of the above analysis should be mentioned. When the 
number of blades is sufficiently large, the coefficients 
C, become functions of Am” * alone (this is always ap- 
for particularly for 
The phase angles then also become 


true compressors, 


proximately 
small k, see Fig. 5). 
functions of A alone for each m (m is only a summation 
index and has nothing to do with the physical situ- 
ation), and one can show that the appropriate extrema 
of each R,,z, over which interval the integral deter- 
mining the corresponding coefficient Q,,2-;(x) is carried 
out, are also functions of A alone, so long as Eqs. (26), 
for example, hold approximately. 

Under these conditions an approximate similarity 
law can be written down, according to which different 
machines operating with the same tip stagger angle and 
resultant tip Mach Number and having the same 
solidity ratio and similar profiles, but having different 
hub ratios and different numbers of blades, will have 


three-dimensional pressure-correction terms related 
through the parameter A. In other words, if 
BD # Bo, yh tH he 
but if 
A, = a a) = hy) — } af a h») = A» = A (56) 
then 
9/3 ie ll => V/1a) 
By “To, iC p oe 6, > OT; Mr, A = 
l1—h 
A l — (7/frr) —_ 
Be %, “Col X08, ; or, Mr, A (D0) 
7 l—h 


This result will clearly be useful in extending the 
range of application of any numerical results obtained. 
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NUMERICAL RESULTS* 


Subsonic Case 


Examples of three-dimensional subsonic _ pressur 
distributions over an axial compressor blade are pr 
sented in Figs. 5 and 6 for two different subsonic ti 
Mach Numbers. 


the pressure distributions according to the ‘‘intuitive 


Also presented for comparison ar 


and the ‘‘modified”’ strip theories. 

In these calculations a blade row with a hub ratio oj 
0.9 and containing 64 parabolic-are blades, each with , 
radial thickness distribution according to Eq. (28 
The axial projection of the chord was 


was assumed. 
taken to be a constant, and under these conditions th 
numbers tabulated in Table 2 can be used directly t 
determine the coefficients Q,,5-;(x). 
profiles the function g(x) [Eqs. (27) and (29) | is given by 


For parabolic-ar 


g(x) = 4U[1 — (2%/cqz) | (58 
consistent with the linearized approximation for thin 
blades. 

In each case a solidity given by L ca, = 2.0 was used. 
The axial Mach Number was held constant at a valu 
of 0.6, with resultant Mach Numbers at the tips of 0. 
and 0.9, respectively. 

In the Figures, values of C,/7o are plotted against 
x C,, at five different radial stations ranging from the 
hub to the tips. The solid lines represent the pressure 
according to the full three-dimensional 
dashed curves correspond to the so-called 
strip theory. The usual, or “‘intuitive, 
strip theory is also plotted at 7/r7 = 0.90 and 0.95 with 
the dot-dash curves. [At the tips (r rz = 1.0) the 
intuitive and modified strip theories are identical. 


coefficients 
theory; the 
‘*modified”’ 


i.e., those 
either 


The three-dimensional correction terms 
three-dimensional effects not included within 
are simply the difference between th 


strip theory 
These differences, 


solid curves and the broken ones. 
or ‘‘corrections,’’ multiplied by B~’’ are plotted in Figs 
7 and 8. The differences that represent corrections 
to the modified strip theory can be correlated approxi 
mately according to the rule given by Eq. (57). In 
other words, according to the theory, these results 
could be applied equally well, for example, to a ma- 
chine having 125 blades and a hub ratio of 0.936 (s 
that A is still equal 1.6), provided the same axial and 
tip Mach Numbers are used, the blades are similar 1! 
shape, and L/c,, = 2.0.3 

It will be noted, from Figs. 7 and 8, that the three 
dimensional effects become relatively more importa 


* The numerical results reported here were largely obtaine 
using the IBM 650 Magnetic Drum Calculator at the Cornel 
Computing Center. 

+ A further example, with the same geometry and Maz 
Mr; = The results are essentially the sam 
in form as those presented here for 1/7 = 0.8 

t Note that application of the similarity rule to a machin 
having significantly fewer blades, e.g., B = 27, would probably 
not be as successful as in the example given, because of the as 
sumption in Eq. (57) that B is quite large 


0.7, was also run. 
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Fic. 5. Subsonic pressure distributions for an axial compressor 
blade row with parabolic-are profiles and B = 64, h = 0.9, L = 
2ez, M = 0.6, Mr = 0.8 at five radial stations ranging from 
hub to tip. Curve A = three-dimensional theory, curve B = 
modified strip theory, curve C = intuitive strip theory T) is 
the thickness ratio at the hub |Eq. (28)], and C, the pressure 
coefficient [Eq. (55 
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Phree-dimensional corrections to the strip theories at 

Curve A = corrections to modified strip theory; 
corrections to intuitive strip theory. Modified and 
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Fic. 6. Subsonic pressure distributions at 7 0.9 and 
Var = 0.6 for a blade row identical to that of Fig. 5. Curves A, 


B, and C have same meaning as in Fig. 5 
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Fic. 9. Supersonic pressure distribution for an axial com 
pressor blade row with double-wedge profiles and B = 64, h = 
0:9, L = Bear, Max = 1.2; Mr = 1.5. Dashed curves give 


Ackeret theory applied according to modified strip theory, solid 
curves are three-dimensional corrections, and their sum gives full 
three-dimensional flow pattern 


as the tip Mach Number approaches unity. This is 
in agreement with the general statement of the theory 
that in the transonic case the flow pattern is inherently 
three-dimensional, to first order in small disturbances. 

Another striking feature of these results is the increase 
in the importance of the three-dimensional corrections 
as one moves in toward the hub. This result is, of 
course, not surprising. The two-dimensional, or strip- 
theory, approximation depends upon the radius being 
large, a condition better satisfied at the tips than at 
the hub. 
two diametrically opposed blades will clearly be stronger 
at the hub. 
sional effects are significant over an appreciable portion 


For example, interference effects between 
In any case, it appears that three-dimen- 


of each blade in the examples treated here. 

The Figures also provide a numerical comparison of 
the relative merits of the intuitive and modified strip 
theories. At these particular subsonic Mach Numbers 
the intuitive strip theory is considerably superior to 
the modified at mid-chord, while the modified strip 
theory is somewhat the better of the two near the 
leading and trailing edges. It is remarkable that the 
intuitive strip theory is so nearly correct at /7 = 0.9 in 
the region of the hub. This result may be due to coin- 
cidence, since the error with respect to the three-dimen- 
sional theory appears to be in the process of changing 
sign at these Mach Numbers. Whatever the reason, 
however, one may conclude that the intuitive strip 
theory can safely be used at higher subsonic Mach 
On the 
other hand, since the modified strip theory is consider- 


Numbers than can the modified strip theory. 


ably simpler to apply, particularly as regards applica- 
tion of experimental cascade results, its accuracy at 
lower subsonic Mach Numbers might be worth in- 
vestigating. 

Since all pressure distributions are symmetric about 


SCIENCES SE 


PTEMBER, 1958 


mid-chord, and since the blades are symmetric, there 


is no wave-energy loss, or wave drag, in the subsonic 
case. This is in agreement with the statements oj 


Rott.° 


Supersonic Case 


Examples of three-dimensional corrections to the 
modified supersonic strip theory for the pressure dis. 
tribution over an axial compressor blade are presented 
in Figs. 9 and 10 for two different supersonic tip Mach 
Numbers. <A blade row with a hub ratio of 0.9 and con- 
taining 64 double-wedge blades, each with a radial thick- 
ness distribution according to Eq. (28), was assumed 
As before, the axial projection of the chord was taken 
to be a constant and the solidity ratio, L c,,, chosen 
to be 2. The axial Mach Number was held constant 
for the two cases at a value of 1.2, with resultant Mach 
Numbers at the tips of 1.5 and 1.8, respectively. 

In the Figures the Ackeret theory for these conditions 
is plotted, according to the modified strip theory, as a 
function of the radial and the chordwise stations. In 
each Figure the ‘“‘positive’’ discontinuity occurring be- 
tween the leading edge and the mid-chord is due to 
the disturbance emanating from the leading edge of 
the blade just in advance of the one in question. The 
“negative” discontinuity at mid-chord in each case is 
caused by the discontinuity in slope of the double- 
The 


arise in similar ways, recalling the “limp body” or non- 


wedge profiles. remaining two discontinuities 
lifting nature of the blade row. 

Applications of the first-order (acoustic) theory in 
supersonic cases where several bodies mutually inter- 
fere suffer severely from the fact that this theory does 
not adequately predict the location in the field of given 
disturbances, although it accurately describes their 
interference problems, 


first-order magnitudes. In 
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Fic. 10. Supersonic pressure distribution for the same blade 


row as in Fig. 9, with Mar = 1.2 and Mr = 1.8. Dashed and 
solid curves have same meaning as in Fig. 9 
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TABLE 3 
Total Wave Drag Coefficients per Blade 


1 D/blade . 43 total torque /blade 
T . ] 2)p U2 792CarrT’ To" * (1 2) p U2 ry VT" 
, a ' Be ss 
Vos 7 » Con _ Cow, = Cow... 2 Ce 
Ti" 2 m i Th" I To" 
1.200 1.500 0.165* 0.170 0.179 0.0895 
1.200 1.800) 0.161* 0.162 0.167 0.1097 


erT a 
* Calculated from f [D. ptt) 4 D(r) \dr 
rH si 


therefore, large errors in location of interference dis- 
turbances occur. This difficulty is often alleviated by 
applying a first-order correction to the location of char- 
acteristic lines (the so-called Lighthill shift). 

In view of this situation, a major part of the differ- 
ence between modified and intuitive strip theories loses 
its significance. The two strip theories differ in two 
respects: first, as to the magnitude of disturbances, 
and second, as to the location of these disturbances in 
the field. The difference in magnitude is obtained by 
replacing the factor 1 wr by | u(r) as the external 
multiplier. The difference in location of disturbances 
in the two theories arises from a difference in effective 
spacing, stagger angle, and slope of Mach lines at every 
radius except the tips. Neither theory, however, 
adequately predicts the location of characteristic lines, 
since each is based on the acoustic approximation. 
Therefore, the difference between the two in this respect 
is moot; a Lighthill-type shift may as well be applied 
to one as to the other. 

In the Figures the three-dimensional corrections to 
the modified strip theory are plotted for five different 
radial stations. Once more, in accordance with what 
we are led to expect from the theory, the three-dimen- 
sional effects are relatively more important at a tip 
Mach Number of 1.5 than at 1.8. Furthermore, their 
importance increases toward the hub. It is interesting 
to note that the modified strip theory in these cases 
actually appears to be superior to the intuitive strip 
theory, particularly with regard to its prediction of the 
drag (see Table 3). 

In the supersonic case, the first-order three-dimen- 
sional effects are largely associated with the swirl of 
Mach ‘‘cones’’ emanating from sources traveling in 
helical paths in a space-fixed system.’ This swirl is the 
more severe, the smaller the radius of the helical path. 
In the present case we are concerned with sheets of such 
sources, and their collective Mach ‘“‘cones’’ form enve- 
lopes that are appreciably distorted, with respect to 
their two-dimensional counterparts, especially near 
the hub. 

The results obtained seem to agree with this inter- 
pretation. Three-dimensional corrections appear to be 
most important in the neighborhood of such discontinu- 
ities, predicted by two-dimensional theory, that arise 
from interference from neighboring blades, and in these 
neighborhoods the magnitude of these corrections in- 
creases toward the hub. Nevertheless, it should be 
pointed out that the magnitude of the three-dimensional 


effects is such as to suggest that they may be largely 
masked by such higher-order phenomena as Prandtl- 
Meyer expansions, shock waves, boundary layers, and 
their respective interactions. In other words, effects 
outside the scope of the present theory may be more 
important, in the supersonic case, than the three 
dimensional corrections here computed. 

The supersonic case, ./,, > 1, presents a consider- 
ably more difficult numerical problem than the sub- 
sonic case, whether the complete formula (15) or the 
analytic expression for the correction terms developed 
in reference | is used. This is because the series con- 
verge relatively slowly when all the X,,, are imaginary. 
The results given in Figs. 9 and 10 were obtained by 
making use of the technique of averaging partial sums, 
and the analytic expressions' for the correction terms 
were used directly. 

The indication of Figs. 9 and 10 that in supersonic 
blade rows three-dimensional effects are probably not 
as important as higher-order effects leads this author 
to believe that the effort required, within the present 
theory, for reliable numerical results will not be found, 
in general, to be justified. Rather, it would seem 
necessary to make use of higher-order theories for 
this case. The essential contribution of the present 
investigation in the supersonic regime is that it provides 
strong justification for carrying out any higher-order 
study of the supersonic blade row within the strip 
theory approximation. 

In the supersonic case, wave-energy losses occur and 
wave drag on the blades appears. Numerical inte 
gration of the pressure distributions given by Figs. 9 
and 10 makes possible the evaluation of the total wave 
drag and total torque coefficients (per blade) according 
to the different theories. The results are listed in Table 
3. It will be noticed that the three-dimensional effects 
in the supersonic case are such as to diminish the total 
wave drag with respect to that computed according to 
The favorable 
It is of 


conventional (intuitive) strip theory. 
effect is stronger at the lower Mach Number. 
some interest that the wave drag according to the modi- 
fied strip theory is closer to the three-dimensional 
value than that given by intuitive strip theory. 


CONCLUSIONS 


A three-dimensional linear theory has been developed 
that can be successfully applied to axial compressor 
blade rows operating in the subsonic and supersonic 
regimes. Although the theory is equally applicable to 
lifting as well as nonlifting blade rows, detailed study 
of only the nonlifting, or thickness, problem has been 
carried out to date. 

Numerical information for determining the eigen 
values of the radial (cylinder) functions under the 
conditions of vanishing radial velocity at the hub and 
shroud |Eqs. (6)] and a large number of blades has been 
provided. The relationship between the full three- 
dimensional solution and the usual strip-theory approxt 
mation, which depends on the radius and the number of 
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blades being large, has been described-—that is, strip 
theory has been shown to correspond to the first term 
of a series development of the flow variables in inverse 
powers of the blade number. Analytic expressions for 


the three-dimensional correction terms have _ been 
worked out,’ and these have been proved to be pro- 


portional to B- (2/5) | 


An alternative strip-theory rule, somewhat simpler 
than the usual one, has been suggested and compared 
numerically with the conventional rule. It is concluded 
that at subsonic Mach Numbers near the transonic 
range the intuitive (conventional) strip theory is su- 
perior to the simpler rule, but in either case three- 


dimensional correction terms are significant. 


Numerical examples of three-dimensional corrections 
are provided in both the subsonic and supersonic cases. 
These corrections can be correlated by a similarity 
rule given by Eq. (57). 


It is further concluded that, in contrast to the sub- 
sonic case, the three-dimensional effects in the super- 
sonic case are probably not as important as higher 
order phenomena, and that higher order studies of 
supersonic blade rows within a strip-theory approxi- 
mation are justified. 


Application of the theory to transonic blade rows will 
be discussed in the second article of this series. 
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Rates by Magnetohydrodynamic Techniques 
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SUMMARY 


heat-transfer 
a mag 
within 


The possibility of reducing stagnation point 
rates on blunt bodies at hypersonic speeds by means of 
The the flow 


netic field is considered modification of 
viscous boundary analyzed. It is 


ind external to the laver is 
primary mechanism which serves to reduce 
alteration of flow external 


‘he existence of a magnetic-fluid bound 


meluded that the 
e heat 
the boundary layer 
layer is indicated for fluids of high electrical conductivity 


transfer is an the inviscid 


iry 
This magnetic boundary layer is seen to parallel viscous bound 
irv-laver theory closely 

» basis of a cursory examination it does not seem prac 


On the 
ticable to reduce significantly the aerodynamic heating load 
m flat plates by magnetic techniques unless the electrical con 


juctivity of air is artificially enhanced 


SYMBOLS 


= velocity gradient parallel to wall 


1 
B = magnetic flux density vector 
E = induced electric field vector 
F = body force vector 
= function related to velocity in Eqs. (A-15) and (B-1 
g = function related to magnetic field, 1. = —g(y 
G = PiL/ Pwhk 
H = magnetic vector 
H = reference magnetic field strength at wall 
h = static enthalpy 
h = heat-transfer coefficient 
= total enthalpy 
| = electric current density vector 
R = thermal conductivity 
= reference length 
p = static pressure 
p = stagnation pressure at stagnation point 
P = Prandtl Number 
Re = Reynolds Number, Re = puol/u 
Ry = magnetic pressure number, Ry = yp. f1?/ pu? 
R, = magnetic Reynolds Number, R,, = opu,to/ 
= velocity component in x direction 
M = reference velocity 
u = velocity vector 
UC = U= yy 
= vertical or ycomponent of velocity 
v, v = Cartesian coordinates along and normal to wall 


} = f. P/ Pe dy 
JI 


viscous boundary-layer thickness 


0 = 

5,,* = magnetic boundary-layer thickness 

6,* = shock-layer thickness at stagnation point 

” = V/}G»,.(X/U,)}"? 

6 = dimensionless total enthalpy ratio, 6 = 
(I — Iw)/(Le — Tw) 

\ = magnetic parameter, \ = ¢,B2,.?/p,.-A 

m = coefficient of viscosity 
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magnetic permeability 


- 7 
1 = kinematic viscosity 
é = R,,'/2 
p = mass density 
a = electrical conductivity 
T = shear stress 
= viscous stress tensor 
¢ = R,,'* f, used in magnetic boundary-layer analysis 
y = stream function, pu = p,,, pv = pu 


Subscripts 


and y components of vector quantities 


= evaluated at outer edge of boundary lave1 


12= x 


relative value, ratio of value with magnetic field to 
value in absence of magnetic field 


= evaluated at wall 


\n asterisk is occasionally used to indicate dimensional quan 


tities; this usage is primarily limited to the discussion of flow 


with large R 


INTRODUCTION 


| acer rivenantemegery concerns itself with the 
mutual interaction of a magnetic field and the 


flow of an electrically conducting fluid. The possi 
bility of using this interaction to advantage in aero 


dynamic problems has been considered in references 

and 2, both of which indicate that significant reduc 
tions in laminar skin friction and heat transfer may be 
achieved under suitable circumstances. It is the pur 
pose of the present report to determine what advantages 
may accrue by using a magnetic field to modify the flow 
in the neighborhood a two-dimensional stagnation 
point, and, in less detail, the flat plate boundary layer. 

The task of evaluating the aerodynamic heating load 
at the stagnation point of a blunt-nosed hypersonic 
vehicle is conveniently divided into two parts —an 
analysis of the viscous laminar boundary layer and a 
calculation of the outer inviscid flow. The high stag 
nation temperatures accompanying flight at extreme 
hypersonic Mach Numbers renders the air sufficiently 
ionized behind the bow shock so that it may be con 
Under these 
field 


as well as 


sidered an electrically conducting fluid. 
circumstances the presence of a magnetic will 
tend to modify both the outer inviscid flow 
the boundary layer. A complete heat-transfer analvsis 
must account for both these effects. 


STAGNATION POINT HEAT TRANSFER 


As indicated in Appendix A, the dependence of the 
heat-transfer coefficient on the outer flow and bound 


ary-layer properties is given by 


h, = Vi 
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WALL 


Fic. 1. Flow about a two-dimensional stagnation point 


In this expression A is the velocity gradient parallel 
to the wall and external to the boundary layer. The 
magnitude of this gradient is determined by the in- 
viscid outer flow. The quantity @,,’, on the other hand, 
is determined by the properties of the boundary layer 
and, in the nonmagnetic case, is independent of A. By 
considering separately the effects of a magnetic field on 
the values of A and @,,’ the overall heat-transfer reduc- 
tion may be obtained. In the course of calculations, 
however, it became apparent that the modification of 
that is, the 


is insensitive to the presence of a magnetic 


the outer flow is the dominant effect 
quantity 8,..’ 
field. 

It may be reasoned that the effect of a magnetic field 
will be such as to reduce the velocity gradient A and 
so too the heat-transfer rate. This is seen by consider- 
ing the equation 


F=jxB 


which gives the body force per unit volume acting on a 
fluid element in terms of the local electric current den- 
sity vector 7 and the magnetic flux density vector B. 
Provided the fluid is electrically conducting, currents 
will be set up in the presence of induced voltage gra- 
dients, that is, gradients caused by the flux cutting 
action of the fluid as it flows past a stationary mag- 
netic field. Assuming a continuous fluid obeying a 
simple Ohm's law relation, the following equations hold: 
j= ck’ and £’=uxB 
so that if only induced voltage gradients exist the body 
force relation becomes 


F = oux8)x8 (2) 


From Eq. (2) it follows that the body force F can only 
have components directed opposite or transverse to 
u; F will always tend to retard the flow. 


THE COMPRESSIBLE VISCOUS BOUNDARY LAYER 


The stagnation-point flow geometry is illustrated in 
Fig. 1, where x and y are taken as the coordinates along 
At some distance 6,* away 


and normal to the wall. 
from the wall, viscous effects cease to be significant 
and, for y > 6,*, the flow is essentially inviscid. At 
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the distance 6,*, which is considerably greater thay 
6,*, the subsonic flow pattern terminates at the boy 
shock. 

In its most general form, the solution of a problem 
in magnetohydrodynamics involves the simultaneoys 
solution of the dynamical equations of fluid motion 
and the Maxwell equations which govern electromag. 
netic phenomena. For the purpose of analyzing th, 
viscous boundary layer, however, several simplifying 
assumptions were made. 

(1) The fluid is electrically neutral, and only induced 
voltage gradients exist. 

(2) The induced magnetic field is negligible com- 


pared to the applied magnetic field throughout the | 


boundary-layer thickness 6,*. 

(5) The fluid, air, is everywhere in thermodynamic 
equilibrium. 

(4) Prandtl Number is constant at 0.71 and the prod- 
uct pu is constant. 
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The first of these is compatible with considering air 
to be a continuous fluid and regarding all solid surfaces 
in contact with the fluid to be at a uniform potential. 
That the induced magnetic field is usually negligible 
will be fully confirmed in the subsequent analysis of the 
and (4+) are felt to be 


outer flow. Assumptions (3) 
justified for the purpose of calculating the relative 
heat-transfer coefficient, /, The variation of elec- 
trical conductivity, o, with temperature was accounted 
for, however, using the theory and electron diffusion 
cross sections of reference 5 for thermally ionized air. 
After suitable mathematical transformations closely 
patterned after reference 4+, the boundary-layer equa- 
tions are reduced to two simultaneous ordinary differ- 
ential equations, as given in Appendix A. These equa- 
tions were numerically integrated for wall enthalpy 
ratios of J, 7, = 1.0 and J, J, = 0, while the stagna- 
tion point conditions were taken to be S,00O0°K. and 
10 atmospheres pressure. These latter values closely 
approximate Mach 25 flight at 100,000 ft. altitude. 
Typical results of these calculations are the velocity 
and dimensionless enthalpy profiles given in Figs. 2 (a) 
and (b) for \ of Oand 1/2. As is seen, the presence of 
an intense magnetic field of approximately S kilogauss 
(which corresponds to \ = 0.5 for Mach 25 and a nose 
radius of 3 ft.) leads to the fuller profiles and thus to 
higher values of skin friction and heat transfer. Prac- 
the alteration of the en- 
summarizes the 


tically speaking, however, 
thalpy protile is negligible. 
results for the boundary-layer calculations giving the 
variation of /”,,, ,and @’,, , with the magnetic parameter 
\. The conclusion then is that, for a given velocity 
gradient A at the outer edge of the stagnation point 
boundary layer, the application of an external mag- 
netic field has practically no effect on the heat-transfer 
rate. Viewed differently, this means that the dominant 
effects of magnetic-fluid interaction in reducing heat 
transfer at a stagnation point occur as a result of modi- 
fying the inviscid flow external to the boundary layer. 


gee 
Fig. 3 


THE OUTER FLOW AT A STAGNATION POINT 


Since the magnetic body force always tends to retard 
the fluid motion, it is to be expected that the velocity 
gradient at the outer edge of the boundary layer will be 
reduced, thus reducing the heat-transfer rate. In ob- 
taining a quantitative estimate of this reduction it is 
necessary to analyze the outer flow, taking into ac- 
count the presence of the magnetic field. To put this 
problem in its proper perspective, it is well to examine 
the significance of the following parameters: 


Re = pl! wu, Reynolds Number 
R, = op, ol, magnetic Reynolds Number 


Ry = p, fy?) pUo?, magnetic pressure number 
6,*, shock-layer thickness 


The equations governing the magnetic and fluid be- 
havior near the stagnation axis (derived in Appendix 
B) for incompressible viscous flow are 
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RADIEt 
Iw/Ie =1.0 
Fic. 3. Summary of boundary-layer solutions 
—(f""’ Re) me fr" en f/2 a 1 = 
Ri(g"? - gg” — Ric*) (3) 
g” = R,,(f’g — g’f) (4) 
where a=xf', HH, = xg’ 
v= —f(y), He = —g(y 


If one assumes the applied magnetic field is vanish- 
ingly small, 1.e., //) ~ 0, then Ry — 0, and the inviscid 
stagnation flow may be approximately represented by 


Fm f. J =y (5) 


By using Eq. (4) it is possible then to trace out the 
behavior of the magnetic lines for vanishing //). In- 
tegrating Eq. (4) subject to the boundary conditions 


gu = —1, gu = 0 


jexp [—(R,, /2)y?]} dy — 
exp [—(K,,./2) y?] (6) 


The magnetic behavior for stagnation point flow is 
then as illustrated in Fig. 4. The significant feature 
is the merging or confluence of the magnetic lines with 
the fluid streamlines.6 In ordinary flow theory it is 
known that viscous effects become negligible whenever 
the vorticity is small. This is the same as saying that 
viscous effects die out as the velocity profile in an or- 
dinary boundary layer merges into the outer flow. 
Defining the viscous velocity boundary-layer thick- 
ness 6,* to be that distance from the wall at which 


uu, = 0.99, gives, for incompressible flow, 
6,*/l = 2.4 Re!’ (7) 


Parallel to this it is reasonable to define a magnetic 
boundary-layer thickness to be that distance from the 
wall at which the vorticity of the magnetic field (which 
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The effect of fluid motion on a weak applied magnetic 
field. 


Fic. 4 


is proportional to the induced current density g”) is 
0.01 of its value at the wall. 
of Ry; = 0 it follows that 


For the preceding case 


6,,.*/l = 3.03/R,,1!/? (S) 


This type of behavior can be expected to persist even 
for Ry, nonzero, in which event magnetic interaction 
with the flow will be largely confined to a magnetic 
boundary layer for large R,,. In present-day practical 
situations, however, F,, will generally be quite small. 
For example, at Mach 25 and 100,000 ft. altitude, o, is 
approximately 4.70 mho/em., and thus R,, is 0.04 
(based on the shock layer thickness 6,* of reference 5, 
a nose radius of one meter, and the velocity m# immedi- 


ately downstream of the shock). From Eq. (8) then 


6.* = 15.1 


which shows that the concept of a magnetic boundary 
layer is invalid unless R,, is of the order of 1,000 or so. 
It follows that for small R,, a good approximation would 
be that the applied magnetic field is virtually undis- 
torted by the electric currents throughout the entire 
shock layer thickness. This simplifies the outer flow 
problem since the values of g and g’ appearing in both 
Eqs. (3) and (4) are known once the applied magnetic 
field is specified. The inviscid outer flow might then 
be solved using Eq. (3) with Re = 
satisfying the appropriate wall and shock boundary 


o; the function / 


conditions. 
by the author. 
static pressure gradient along the wall is unaffected by 


Such an analysis has not been carried out 
In reference 2 it is assumed that the 


the presence of a magnetic field. With this assumption 
the reduction in A may be directly computed by evalu- 
ating Eq. (A-9) outside the boundary layer and for 
i= 
about 7 per cent. 


0.5 gives an overall heat-transfer reduction of 


Although there seems to be no immediate aerody- 
namic need to consider the case of large R,,, a short dis- 
cussion of large R,, phenomenon is called for as an aid 
situations 


to general understanding. Besides, 


may 
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well occur in which it is no longer valid to completely 

disregard the distortion of the magnetic field. 
The appropriate stagnation point equations to bg 

solved are Eqs. (3) and (4). It is convenient to intro. 

duce the simple change of variables 


= K,,''*y, o = K,,"*} 


which gives 


[( —R,,) Rely’”’ ge’ re ¢”? —— 
Rui Rn(g”? — gg” ( G 
o” = gg — g’¢ (10 


The wall boundary conditions are | 
/ 
a= = iL, Lu _ 0, Su = Su = 


while the outer conditions are 


9 


~ ’ , e 
E> ono, go ~ 1.0, g °c (11 


the requirement that g’? ~ c? as > © expresses the 


confluence of magnetic and fluid streamlines. <A solu 
tion of Eqs. (9) and (10) was obtained by expressing 


g and g as series in powers of RiR,,; i.e., 


yg = $0 + RnR, + (RiRy,)*¢2 ; aa f 


o = + RoR, + (RoR) +... § 


The resulting solution for the velocity profile correct to 
the first power in R7/R,, is shown in Fig. 5 for RR, 

0.345 and R,, Re = 1.0. For this value of RR, the 
wall shear is reduced 50 per cent while the heat-transfer 
coefficient is reduced 15 per cent for Pr = 0.71. As 
R,, decreases the magnetic layer thickness becomes ma 
terially greater than the viscous boundary layer. It 
is conceivable that the circumstance may arise in which 


6 —.. <= 6. 


In this event, the conventional viscous boundary layer 











7 With this boundary condition it follows that RyRy = 
oB;?/pA as RyRm = op.?He?l/ pio = oc By2l/pug and ue* = Ax* or 
te = (Al/uy)x; consequently, //u) = 1/A 
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Fic. 5. Constant property stagnation point velocity profile for 


Rn = R.. 
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in which inertia and viscous forces are comparable) 
would be buried beneath a magnetic layer (in which 
magnetic and inertia forces are comparable) which in 
turn lies beneath the outer near-potential flow region 
where inertia forces predominate). 

One other interesting consequence of large R,, is 
that a transverse pressure gradient exists across the 
magnetic boundary layer. This may be seen by con- 
sidering Eq. (B-6) when written in the £ coordinate 
system—namely, 


p2 = (RR, 2) (g’? — c“-) + (1/2) (133) 


Thus, Ape = po, — Po, w = (RuRn/2)c? (14) 


From Eq. (B-1) the static pressure p is given by 
Pp = po — pily) — poly)x* (15) 


and since 
Api Pi.e — Pi = (7.7/2) =f, /) = 
( ~0 for large R,, 


then Ap = —Apox? = —(RyR,,/2) cx? (16) 


According to Eq. (16), then, the static pressure at the 
wall is greater than the pressure at the outer edge of 


the magnetic layer. 


FLAT PLATE HEAT-TRANSFER RATES 


In the analysis of compressible stagnation point heat 
transfer it was noted that the magnetic effects on the 
outer flow were dominant. For a flat plate at hyper- 
sonic speeds the situation is reversed since outside the 
viscous boundary layer the air is cold and not ionized. 
Consequently, the electrical conductivity is zero out- 
side the boundary layer and no magnetic-fluid inter- 
action occurs there. Since the degree of ionization ts 
primarily a function of temperature, the magnetic 
interaction with the fluid in the boundary layer will be 
greatest where the static temperature maximizes. By 
examining the magnetic parameter R,/R,, it is possible 
to obtain a quick estimate of how significant the mag- 
netic forces are. If the wall temperature is limited to 
some reasonable value, say 1,600°K., then the maximum 
static temperature in the boundary layer will be ap- 
proximately 2,750°K. according to reference Sf for a 
free-stream Mach Number of 25. 

Evaluating p, o, and mw at this point for a pressure 
corresponding to 100,000 ft. altitude gives 


RnR, = cB2?l/puo = 0.002 


ior a field of 10 kilogauss and an / of one meter. 

This discouragingly small value for R,R,, is largely 
due to the relatively low peak static temperature in the 
boundary layer and thus electrical conductivities of 
only | 10,000 as great as at stagnation point conditions. 
It seems reasonable to conclude that there is little ad- 
vantage in magnetohydrodynamic methods in reducing 

| Although based on perfect gas assumptions, the method of 


this reference should be sufficiently accurate for order of magni 


tude purposes 
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flat plate heat-transfer rates. If practical means of 
artificially increasing the electrical conductivity of 
the boundary layer are found, then this conclusion 
could well be reversed. 

APPENDIX A— THE COMPRESSIBLE VISCOUS 
STAGNATION POINT BOUNDARY-LAYER EQUATIONS 


According to reference 7 the magnetohydrodynamic 
equations for the steady flow of a continuous electri 
cally conducting fluid with zero net charge are as 
follows. 

Fluid equations: 


Continuity V-(pu) = 0 A-1 
Momentum p(Du Dt) = -Vp+V-74+F (A-2) 
Electromagnetic body force F = jxB (A-3 


Energy p(DI/Dt) = V-(u7) + V-(RVT) 4 


DiGge uxB) (A-4 
Electromagnetic equations: 
VB=0, j= Vx (A-5 
and the constitutive relations 
B =z 
. (A-6 


J = a(E + E’) = cE T ou x BI 


If it is assumed that no applied or polarization volt 
ages exist that is, E = 0—then the electrical terms dis 
appear from the energy equation. This then corre 
sponds to the case where no energy is being added or 
extracted from the fluid by electrical means. The 
magnetic effects are confined to retarding the flow and 
dissipating energy of motion into heat. 

In general, the electrical currents flowing in the fluid 
will give rise to an induced magnetic field which will 
distort the applied magnetic field, that field which 
would exist if the fluid were an electrical insulator. 
However, since the viscous boundary layer is thin and 
thermally ionized air is at best a poor electrical con 
ductor it is permissible to neglect the induced field 
compared to the applied field. This permits the simphi- 
fication B,; = 0, B., = constant, and the body force be 


comes 


F, = —oB.u, F, = 0 (A-7) 


Thus the only modification to the conventional bound- 
ary-layer equations is a —oB*u term in the momentum 
equation for the direction parallel to the wall. 

The required boundary-layer equations for the stag 


nation point then are 


Continuity {O(pu) Ox] + [O(pv) Oy] = O (A-S 
Momentum equation along wall 
pu(Ou Ox) + pv(Ou Oy) = 
—(Op Ox) + (O Oy) [u(Ou Oy) ] oB*u (A-9 
Op oy = 0 A-10) 
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pu(Ol/ox) + pu(Ol Oy) = (0 Oy) [(u/Pr) (OF/Oy)] 


(A-11) 
Introducing the stream function 
Pu Po = Wy, pl! Pu —yr 
and the transformation of coordinates 
Xanax, ¥= | Bicls 
gives 
UUx + VUy = (oc/p)UUey + 
(By?/p) (o,U. — cU) + wy (GUy)y (A-12) 
where 
G = pt/puky = constant, U = yy, V = —yx 


Writing the equation of state in terms of enthalpy as 
p = pzRh = pzRI 
where z is a modified compressibility factor leaves 
UUx + VUy = UUex(l/Te) (2/2) + 
(B.?/p) (o.U, — cU) + »,(GUy)y 
The transformation 
7 = Y/(Gv.(X/U,)]!/? 
y = U,.|Gr,(X/U,)]'" fin) 
and assuming 
yields 
fr 4g” =f? — (1/I,) (2/%) — 
(B.?/A) [(o. — of’) p| 
where p, o, 2 2, and J J, are to be functions of n only. 
The energy equation transformed into the X ¥ plane 
is 
Uly + Vly = ww [(G/Pr)Ty]y 
(I — Iy)/ULe — Iw) = O(n) 


and writing 
gives 6” + Prfo’ = 0 

For convenience the magnetic term in the momentum 
equation is written as 


(Bx?/A) (oe pe) \(pe/p) — [(o/ ae) (p p.) If} 
or A[(Z/I,.) (2/2) — Jf] 
where \ = B.*c./Ap, (A-15) 


The parameter \ is then a measure of the importance 
of the magnetic interaction term as compared to the 
inertia or viscous terms. The two ordinary differ- 
ential equations for f and @ are to be solved simultane- 


ously subject to the boundary conditions 


»=0: f=f’ =0,06=0 
f’ > 1.0, 6 => 1.0 


Once these equations are numerically integrated the 
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skin-friction and heat-transfer coeflicients are obtained 


from the relations 


Cy = 2re/ putts? = (2f5"/x) (¥e/A)? (A-16 


hy = [Rw(Ol Ov) w/e — Le)] = Ru(A/vw)'!? Oy’ (A-17 


APPENDIX B—Exact EQUATIONS FOR’ CONSTAN} 
PROPERTY MAGNETOHYDRODYNAMIC STAGNATION Porn] 


FLow 


The appropriate basic equations are (using an aster. 
isk to denote dimensional quantities) 

V-u* =0, V-H* =0 

pu*-V*u* = —V*p* + gigs” x H*) + uwV**u* 

J* = ou,(u*x H*) (V*x H*) 


Substituting the dimensionless quantities 


x=x*/l, y= y*/l, w= u*/ue 
p = p* ply”, H —_ Df TI, 
gives 
u-Vu = —Vp+ Ry (Vx) xH+ (1 Re)V7u 
Vii = Kk, @ =f) 
Vu =V-H=0 


Assuming the fluid and magnetic stream functions to 
have the form 


vr = xf(v), Yu = xg(y) 
it follows that 
u=xf', HH, = xg’ ] 
v= -—-f, HM, = -g (B-1 
p = po — pily) - po(y)x°f) 
The previous set of equations then becomes 
((—f'"") /Re] — ff” +f’? = —Ruge” + 2p. (B-2 
ff’ = pi’ — (f"/Re) (B-3 
po’ = Rug’g” (B-4 
g” = R,(f'g — g’f) (B-5 
Integrating Eq. (B-4) gives 
pe = (Rug’?/2) — (RuR,,/2) c2 + (1/2) (B-6 


where —RrR., (¢? 
tion. 


2) + (1/2) is a constant of integra- 
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Tt These functional forms coincide with those in reference 2. 
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Some Effects of Curvature on Frames 


ROBERT W. WESTRUP* anp PHILIP SILVER** 


North American Aviation, Ine. 


SUMMARY 


Induced radial loadings due to curvature can cause drastic 


changes in the stress distribution on the cross section of aircraft 
frames from that predicted by elementary bending theory 

In the present paper, the exact theoretical solution for the 
stress distribution is reviewed for simple plate elements. An 
ipproximate solution for the stress distribution which is appli 


Gor rd 


between the two theories is indicated for the basic 


cable to more general cross section shapes is then developed 
igreement 
model where they can be compared 
js defined, and design charts showing its variation with the cur 


An “effective area’’ concept 
vature parameter are presented 

Experimental investigation of the stress distribution and ulti- 
mate crippling capacity of typical frame specimens is outlined 
Excellent agreement was obtained between test results and the 


ipproximate theoretical solution 


SYMBOLS 


Parame.ers 


A area 
h width 
C arbitrary constant of integration 
E = modulus of elasticity 
h total frame depth 
m = virtual bending moment 
V = bending moment 
r = load 
R = radius 

= coordinate along web 
t = thickness 
V shear 

running load 

\ coordinate along flange 
J = radial deflection 
a, nondimensional curvature parameters 
B = curvature parameter 
6 = reference deflection 
€ = unit strain 
Mm = Poisson’s ratio (assumed equal to 0.3 
” = efficiency 
o = reference rotation 
a = circumferential stress 


Subscripts 


F = flange element 
i = lip element 
i = web element 


INTRODUCTION 


Pparapene TESTS conducted on missile tank frames 
indicated failing loads considerably lower than the 
predicted values. Subsequent testing with more ex- 
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* Design Specialist, Structures Group, Missile Development 
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tensive instrumentation showed that the variation was 
due to nonuniform stress distribution, resulting in loss 
of effectiveness of outstanding flange elements of the 
frame caps. The general problem has been considered 
by Timoshenko! and Seely,” and solutions for curved 
beams of simple cross section are given in references | 
and 2. However, these results have somewhat limited 
application. In the present paper, an approximate 
solution for the stress distribution in curved beams of 
more general cross-section shape is presented, together 
with the results of an experimental investigation to 
check its validity. The exact solution for stress dis- 
tribution on simple plate elements under various edge 
support conditions is also reviewed for design use and 


for comparison to the approximate theory. 


NATURE OF STRESS DISTRIBUTION 


The general characteristics of the stress distribution 
on a curved frame cap with outstanding flange ele- 
ments is illustrated in Fig. 1. It is seen that the stress 
is highest at the vertical web element and decays along 
the flange element with distance from the web line. 
This decay results from distortion of the flange cross 
section. The circumferential produces an 
equivalent radial loading which in turn causes transverse 
A segment 


stress oa 


bending deflections of the flange elements. 


of the frame cap, as depicted in Fig. 1, will shorten 
some amount under the indicated compressive load- 
However, fibers out along the flange element can 


(1) by 


ing. 
match this general deformation in two ways: 
actual fiber strain, assumed to follow elastic behavior; 
and (2) by deflecting to a different radius of curvature, 
thus matching the general deformation without strain- 
ing. 

Therefore, the transverse bending deflections vary 
the radius of curvature along the flange element and 
result in a nonuniform circumferential stress distribu- 


tion. Previous investigations!: * have defined the rela- 
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Fic. 1. Curvature effects 
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tionship between these two quantities as 
or = ow — E(y/R) (1) 


It is still necessary to relate the flange distortion to 
the circumferential stress distribution so that Eq. (1) 
may be evaluated. If the restraint along the length of 
the frame cap is neglected,* the deflection can be evalu 
ated as the elastic curve of a transverse strip. 

The induced radial loading on a transverse strip of 
unit width, illustrated in Fig. 1, can be expressed as 


w= o,t/R = (owt/R) — (Ety/R?) (2) 
EXxAcT SOLUTION 


Stress Distribution 

For a frame cap where the flange can be represented 
as a simple plate element, the classic equations for bend- 
ing curvature of a plate and relationship of running 
load to moment may be used to modify Eq. (2) to an 
expression containing only y as the dependent variable. 


d*y/d,* = M/D 


where DP = £F/I20 — p*) 


w= d°M/dx? = D(d‘y/dx*) 


Substitution of these relationships into Eq. (2) and re- 
arranging yields 


(d+y/dx*) + kty = k*(Row/E) (3) 


where kt = 12(1 — p?)/tp?R? 


The general solution to this differential equation may 
be written in the form 


y = —(Row/E) [(C, sinh Bx + C. cosh Bx) sin Bx + 
(C; sinh Bx + Cy cosh Bx) cos Bx — 1] (4) 

where B4 = k4/4 = 3(1 — p?)/t?R’ 
Substitution into Eq. (1) yields the general equation 
for circumferential stress distribution 

* This is perfectly permissible if the loading does not change 
along the length of the frame cap. It should be an acceptable 
approximation for most frames with normal rates of loading vari- 


ation. 


to structural design or analysis. 
area,” 


the flange element.t 
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(C; sinh Bx + Cy cosh Bx) sin Bx + 
(C; sinh Bx + C; cosh Bx) cos 34 5) | 


o OV 


The arbitrary coefficients must be evaluated by con 
sidering the appropriate boundary conditions. Mos 
frame caps would be idealized as cantilever plates 
however, the solutions for simply supported and fixed 
edges will be included for possible application to box | 
frame or wide curved panel configurations. Thy 
boundary conditions and resulting constants are listed 


as follows: 


Cantilever plate 


x 0, y = (- C; | 

x 0, dy/dx = 0; Cs —C; 
d*y 

, ; ax* °% " 


= sinh? Bh) + sin? Bd 
cosh? 6b + cos? Bb 
dy 7 
x = 5b, de = 0; C; 
7 sinh Bb cosh Bb + sin Bb cos Bb 


cosh? Bb + cos? Bh 


Sim ply supported edges 


x=0, y = 0); C; = 1 

vr = 0, 2 0; C, = 0 
dx* 

\ b ¥y QO; Cy 


(cosh Bb — cos Bb) sin Bb o 


sinh? Bb + sin? Bb 


(cosh Bb — cos Bb) sinh Bb 


sinh? Bb + sin? Bb 


x =0, y = (): C; = 1 
v= 0, dy/dx = 0; om «€; 
x=b y = 0: Cc, = 


sin 8b — sinh Bb 

sinh 6b + sin Bb 
dy/dx = 0; C3; = 

cos 8b — cosh Bh 


sinh Bb + sin Bb 


Effective Area 


Eq. (5), defining the circumferential stress distribu- 


tion, is of theoretical interest, but is not easy to apply 


A more convenient 


way of expressing the result is to define an ‘‘effective 


such that the product of this area and the maxi- 


mum flange stress ow gives the correct load carried by 


The effective area may be de- 


fined by an “‘efficiency factor’ y, such that 


t This is completely analogous to the ‘effective width’’ con- 


cept applied to compression panels with buckled skin 
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Eg. (5) may be used with Eq. (9) to express the effi- 


ciency as 


n= [1 V3 — w2)] (VWER/D)C, 10 


where C; is a function of the previous constants (C, 
C2, C3, Cs) 
gration of Eq. (9). 

This constant has been evaluated for the three flange 


and additional terms arising from the inte- 


support conditions as follows: 


Cantilever plate 
sinh 8b cosh 6b + sin 8b cos Bb 


cosh? Bb + cos? Bb 


Simply supported edges 
(sinh Bb + sin Bb) (cosh Bb — cos Bd) 


sinh? Bb + sin? Bh 


Fixe (l edges 
2(cosh Bb — cos Bb) 
sinh 6b + sin Bb 


Efficiency curves for the three support conditions are 
plotted in Fig. 2 against the frame parameter (RK /}’. 


APPROXIMATE SOLUTION 


The exact solution to the differential equation de- 
scribing flange behavior may be found quite readily for 
the case of simple plate elements. However, many 
practical frame cross sections will have perturbations 
such as bulb or lip elements at the tip of the flange, or 
will be built up in such a manner that one may no 
longer assume a condition of fixity at the web center- 
line. In such cases, the general solution to the differ- 
ential equation would be difficult, if not impossible, to 
obtain; therefore, an approximate solution has been 
developed which is believed to give good accuracy in 
the practical design range. A general deflected shape 
is assumed, and the stress distribution and flange effi- 


ciency are evaluated on this basis. 


Stress Distribution 


A typical frame cap of rather general nature 1s illus- 
trated in Fig. 3(a). A bulb or lip area .1;, is assumed 
at the tip of flange, the effective point of fixity is as- 
sumed some distance fy down along the vertical web 
element, and the thickness of flange and web elements 
are considered to be different. 

The induced radial loading on a transverse strip of 
unit width is indicated in Fig. 3(b). The lip element 


provides a shear force 


P, = 0,A,/R (11) 


ATURE ON FRAMES bu 


and the running load over flange element is expressed as 
Wp ot R 12 

It is assumed that the elastic curve of the flange ele 
ment is a portion of a sine wave, with 6 as the deflection 
at tip, and that bending of the web element causes a 


rotation ¢ at root of flange. 
The equation for total deflection then becomes 


y 6/1 — sin (mx 2b, + ob; x 13 
where «is now measured from the tip of flange. 


With this assumed deflection curve, the general stress 


distribution, Eq. (1), becomes 
o Cv -— (E R) )6[1 a sin ( 7x 2b; T 


dolby — x){ l4 


and the induced radial loading may be evaluated as 


P, lon = (i R) (6 + ob; | (.] R lo 
Wy (tre R) (on — (E/R) )6/1 _ 
sin (rx/2br)] + (bp - x)}) 16 


The transverse bending moment distribution may be 
obtained by successive integration of the loading equa 


tions 


V; Pi + wredx, Vy 0 


My | V pdx, Mnry AM, : bE 17) 


with the result 


AM; (Etrd GR*)x? + 
l(owtr 2R) = (fet; 2h?) (6 7 dbp) x T 


[(owAz/R) — (2£.A,/R?*) (6 + ober) 4 
(2tvEébp/ rR?) |x — (AE6bp*tp / r?R?*) sin 
(wx/2br) (18 


The reference deflections 6 and ¢ can be evaluated by 
a number of techniques; the theory of virtual work is 
probably one of the most direct and will be used here. 
The virtual bending moments resulting from a_ unit 
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load applied at the tip of the flange and a unit mo- 
ment applied at the junction of flange and web elements 


are, resp ctively, 
my x Me = | 


and the reference deflections are 


p 12(1 — p?*) i’ 
0 : 
Etp® J0 
12(1 — yp’) 
Etw? ‘ 


Mpmpdx 


bd 71 
| My my ds 


These integrals may be evaluated and put in the 


following nondimensional form: 


6 b; 
(on kk) (R bp) [1 


ri) [(1/2) + (41 
(ow E) (Ribp) [I 


Eqs. (19) and (20) could be solved together to ex- 
tract Y and Y independently, but the resulting general 
expressions would be very cumbersome. If numerical 
values are inserted at this point, the solution is greatly 
simplified. 

The general stress equation can be rewritten in terms 
of the nondimensional parameters Y and Y as 


by) | = 


Oo, on )-(X 
Y{1 — (x/br)] (21) 


1 — X [1 — sin (7/2 


The relative stress distributions for the frame speci- 
mens tested in this investigation have been calculated 
and are plotted in Fig. 8. Strain gage points are in- 
cluded for comparison. In these calculations, the ef- 
fective point of fixity was taken as the centerline of the 


rivets through the web elements of cap angles. 


Effective Area 
An “effective area”’ and “efficiency factor’’ can be de- 
fined in the same manner as before for both the flange 


and lip elements. 


bd J 


Flange element nr = | otpdx/owtrby 


0 


Lip element nL = O1/ow 


With the help of Eq. (21), these expressions may be 


written as 
ne = 1 — [1 — (2/m) JX — (1/2) Y 
m=1—-X —Y 


The efficiency for a simple cantilever flange, as deter- 
mined by the approximate solution, is plotted as a 
dotted line in Fig. 2, where it may be compared with the 
exact solution. It is seen that for values of tk}? in 
excess of 0.5, the agreement is very good. The practi- 
cal design range would normally be above this point, so 
the approximate method is considered to be quite 


satisfactory. 


RVATURE 


(1 Ss) + (] 3) ° 


12(1 — w?)] (tryR/bp?)? + (1/3)- (C1, 
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Ay — [(11/120) + (1/3)-C1,/-1¢)]) 
: 19 
Apr) + [(1/8) + (16/24) — (2/37 
Ap)|] — (1 2) — (2/r) + (4 ry) + (1 1; 4 2») 


Design Charts and Application 


The resulting equations for stress distribution or 
efficiency are quite cumbersome to apply to routine 
design or analysis work, so a series of design charts 
have been prepared to expedite the application of the 
theory Some sample design charts for the efficiency of 
flange and lip elements are presented in Figs. 4 and 5. 
Note that the efficiency of the lip element may become 
negative for small values of frame parameter. 

The manner in which various practical frame cross 
sections may be idealized to the general configuration 
treated in this analysis is indicated in Fig. 6. The un 
symmetrical channel section [Fig. 6(c)| involves bend 
ing of the web. Test results indicated an effective 
length of web element equal to 0.4 of the total depth to 
be a reasonable empirical value. 


ADDITIONAL EFFECTS 


Crippling 

The stress distribution predicted in the previous 
section is based on the assumption of elastic behavior. 
If this same general trend persists well into the plastic 
range, it can seriously affect the ultimate crippling load 
capacity of the frame cap. A crippling test was per- 
formed on one of the frame specimens, as described in 
investigation. 


the discussion on the experimental 


While the relative strain distribution was somewhat 
different than in the elastic range, a large reduction in 
The 


actual ultimate load was about 30 per cent less than 


load capacity of the flange element was noted. 


what would be predicted for a straight flange of the 
same section, assuming a crippling cutoff equal to the 
yield stress of 40,000 psi. An ultimate load predicted 
by using the same crippling cutoff with the effective 


flange area is within 2 per cent of the measured value. 





572 JOURNAL OF THE AERO 


Lateral Instability 

Frame curvature reduces generally the effectiveness 
of flange elements, with the region near the free edge 
being most severely affected. This means that the 
effective lateral bending stiffness is reduced, probably to 
a greater extent than the axial stiffness. This should 
certainly be considered in predicting a critical stress 
in the lateral instability mode. No rigorous theo 
retical treatment is attempted here; however, approxi 
mate corrections have given fair agreement with ex- 
perimental data. The approximate correction as 
sumed the effective flange thickness distributed in the 
same manner as the relative axial stress distribution; the 
effective lateral moment of inertia was calculated on this 


basis. 
EXPERIMENTAL INVESTIGATION 


Circumferential Stress Distribution 


Four curved frame specimens were used to investi- 
gate the accuracy of the stress distribution equations. 
Three of these specimens were of identical cross section 
and formed to various curvatures (inner cap radius of 
15, 30, and 
duction frame with lipped flanges and an inner cap 


15 in.). The fourth specimen was a pro- 


radius of 25 in. Axial strain gages were installed at 
about the midsection of each specimen to completely 
define the normal stress distribution over the cross sec- 
tion. The cross section of the specimens and the ap- 
proximate strain gage locations are illustrated in Fig. 7. 

Each specimen consisted of a pair of identical frame 
The 


specimens were tested in a special loading jig which ap- 


segments* separated by a skin panel 10 in. wide. 


plied a pure moment to each end of the specimen, such 
that the free flanges were in compression. 

Strain readings were taken at several load incre- 
iments so that a linear load-strain relationship could be 
established in the elastic range. The relative strains 
at gage points on the flange elements are plotted in 

* The circumferential length of the specimens was about 48 in., 
so the included angle varied from about 60° to 180°, depending 
on the curvature. 


SPACE 


SEPTEMBER, 1958 


SCIENCES 
Fig. 8. These data are the average of the two gage 
readings at corresponding points of left and _ right 
flanges, divided by the maximum strain at the web 
line. The strain gage data may be compared with the 
theoretical stress distribution; excellent agreement js 


indicated. 


Ultimate Load Capacity 


After the previous investigation was completed, the 
frame specimen with greatest curvature (Kk 15) was 
loaded until ultimate crippling failure occurred. 
Strain readings were taken at regular load increments 
the last set of data was obtained at 9S per cent of ulti- 
mate load. The strain was extrapolated to ultimate 
level and converted to stress distribution by use of a 
typical stress-strain curve for the material. The stress 
distribution is plotted in Fig. 9. The ultimate load 
capacity of the flange was obtained by numerical in- 


tegration of this stress distribution. 


CONCLUSIONS 


The effects of curvature can be very important in 
aircraft frame design and should be considered in pre- 
dicting applied or allowable stresses. 

An approximate method is developed to determine 
the applied stress distribution on frame caps of general 
cross section. It shows good agreement with test re- 
sults and with the exact theoretical solution where com- 
parison is possible. 

Efficiency charts may be constructed from the ap 
proximate solution to expedite detail analysis and assist 
in the selection of efficient cross-section proportions. 
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- Thermal Buckling of Solid Wings of Arbitrary 


Aspect Ratio’ 


JOSEF SINGER* 
Polytechnic Institute of Brooklyn 


SUMMARY 


The thermal buckling of solid wings with all edges free is 
First wings of large 
ispect effect 


double series of orthogonal polynomials is assumed for the de 


the Ravleigh-Ritz method 
the 


nalyzed by 
ratios are treated and end is neglected. A 
flection function, and by the use of successively increasing num 
bers of terms the predominance of the simple torsional distortion 
is established for a number of typical temperature distributions 
for wings of constant thickness and of parabolic cross section 
The effect of variation of aspect ratio in this analysis is also ex 
umined. The validity of the conclusion about the predominance 
of the simple torsional distortion is then verified by a buckling 
nalysis for a wing of constant thickness and small aspect ratio, 
using an approximate solution of the thermal-stress problem of a 
finite rectangular plate 

Having been established as a good approximation, the tor 
sional analysis is now modified to include wings of small aspect 
ratio also. It is shown that the effects of chordwise normal 
stresses and shear stresses cancel out, and that the critical tem 
perature of a finite wing can be obtained from that of an infinite 
wing by multiplication by an ‘‘end effect factor.’’ This factor is 
calculated for various temperature distributions for wings of 
constant thickness, parabolic cross section and double wedge 
shape, and shown to be approximately the same for different 
typical temperature distributions and nearly unaffected by the 
shape ol the cross section. The end effect is found to be signifi 
cant (more than 10 per cent for aspect ratios below 6), but can be 


included with very little additional work 


SYMBOLS 
semispan of wing, ft. or in 
= coefficient of a typical term of the double series 
expansion of the general deflection function 


semichord of wing, ft. or in 


( = effective polar moment of the cross section, in 
( = constant defined by Eq. (23 

da, Le, = constants 

D = bending stiffness = Fh*/12(1 — v?), Ib. in 

E = Young’s modulus, psi 


= function of y only 
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F function of y only 
FE = constant defined by Eq. (24 
function of x only 
G = shear modulus psi, or arbitrary function of x 
h = thickness of wing, ft. or in 
h maximum thickness of wing, ft. or in 
H constant defined by Eq. (25 
twist per unit length, rad. per in 
i, constant defined by Eq. (26 
Ni = constant defined by Eq. (27 
stress resultant in x direction, lb. per in 
N stress resultant in y direction, Ib. per in 
r = Legendre polynomial of degree n 
Q," ultraspherical polynomial of order r and degree n 
\ shear stress resultant, Ib. per in 
7 temperature of wing, “R 
‘i? adiabatic wall temperature, °R 
1 initial temperature, °R 
/ temperature parameter, °F 
l strain energy, Ib. in 
J potential of external forces, Ib. in 
deflection, in 
a coefficient of thermal expansion, ft. per ft., °F 
r end effect factor defined by Eq. (34 
n nondimensional coordinate 
v Poisson's ratio 
a stress in x direction, psi 
o = stress in x direction in wing of infinite aspect 
ratio, psi 
a stress in y direction, psi 
Try = shear stress, psi 
V = stress resultant function 
é nondimensional coordinate v/a 
INTRODUCTION 
AERODYNAMIC HEATING of solid wings or tail 


TT. ; 
surfaces, as frequently used in missiles, causes 


transient chordwise temperature variations which give 
rise to spanwise thermal stresses. When the tempera 
ture of the leading and trailing edges is higher than that 
of the center portion of the wing, as in accelerated 
flight, buckling may occur. 

loss of 


the 
torsional stiffness caused by aerodynamic heating in 


Budiansky and Mayers! investigated 
solid wings of symmetrical double-wedge cross sections 
accelerated to supersonic speeds, and found that the 
losses may be severe enough to produce torsional buck 
ling. The effect was also demonstrated experimentally 
by Vosteen and Fuller.” 

With the assumption that the displacement of the 
wing, when it buckles under thermal stresses, is a simple 
torsional distortion, and neglecting end effects, Budian 
sky and Mayers derived a simple formula for the effec 
When the ef 
fective torsional stiffness vanishes, the wing buckles 


tive torsional stiffness of a solid wing. 
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Fic. 1. Transient temperature distribution for a wing of 


parabolic cross section (¢/c ratio 0.02) accelerated instantaneously 
toa Mach Number of 3 in the stratosphere. 


and then the formula yields the critical stress. Essen- 
tially the same formula was also derived simultaneously 
by Hoff* and Bisplinghoff' by different approaches. 

The torsional analysis yields a convenient expres- 
sion for the critical temperature parameter (which in- 
dicates the onset of buckling) for a wing subjected to 
thermal stresses ¢, = ak1\F(y): 


eh 
Ter = —GC ak: | hy* F(y) dy (1) 
/ b 


The critical temperature parameter calculated from 
Eq. (1), corresponds to the particular temperature dis- 
tribution which induced the stress distribution F(y). 
The magnitude of this temperature distribution is de- 
fined by a temperature parameter 7). Buckling will 
take place if 7)., for this distribution is equal or less 
than the actual 7}. 


On heating, or cooling, the temperature distribution 
of the wing varies with time both in magnitude and 
shape, as can be seen, for example, in Fig. 1, which 
shows the temperature distribution in a solid wing 
of parabolic cross section accelerated instantaneously 
to a Mach Number of 3 in the stratosphere (see Ap- 
pendix B of reference 5 for calculation of the tempera- 
ture distribution). The whole time range has therefore 
to be investigated for possible buckling by repeated ap- 
plication of Eq. (1). 

However, the reliability of the analysis yielding Eq. 
(1), depends on the validity of the assumption that the 
lowest buckling mode is a simple torsional distortion, 
or that the buckled shape can at least be approximated 
closely by a torsional deformation. If this assumption 
is valid, a general buckling analysis, without any re- 
strictive assumption on the displacement function, 
should demonstrate the predominance of the term cor- 
responding to torsional deformation. Such an analysis 
is carried out here. 

The accuracy of Eq. (1) depends also on the magni- 
tude of the error introduced when the effect of the free 
edges is neglected. An approximate investigation of 
the free end effect is carried out in reference 1; in the 
calculation the chordwise normal stresses are ignored. 
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A more accurate analysis, which includes the chordwise 
normal stresses, is carried out here and it is found that | 
the effects of shearing stresses and chordwise norma! 
stresses cancel each other out in the case of torsional 
distortion. 

In the present calculation a wing with all edges fre 
is considered. This idealized model is assumed to rep. 
resent a typical solid wing as used in missile construc- 
tion. The free-free boundary conditions are preferred 
to those of a cantilever (which imply rigid supports), 
since here the supports do not provide much restraint 
to warping or change in shape of the wing. 


BUCKLING ANALYSIS FOR WINGS OF LARGE ASPECT 


RATIO 


In order to simplify the general buckling analysis 
which is to demonstrate the predominance of the tor- 
sional deformation, wings of large aspect ratio are con- 
sidered first. 
at the tips can be disregarded and the thermal stresses 


For such wings the local stress effects 


can be approximated by those of a wing of infinite as- 
pect ratio. 

The wing to be analyzed is shown in Fig. 2, and here 
the semispan a is very large compared to the semichord 
b. Due to aerodynamic heating the wing is subjected 
to a chordwise temperature variation 7) = 7(y), which 
can also be written as 7) + 7)f(y). The resulting 
stress distribution is (see Section 1.8 of reference 6) 


or = akT\F(y), a, = 0, ty = 0 (2 


Now if the simple torsional distortion were an exact 
buckling mode, it would obey the differential equation 
of the rectangular plate subject to the thermal stress 
distribution of Eq. (2) and satisfy the free-free bound- 
ary conditions. The deflection function of a simple tor- 
sional distortion is 


w= kxy (3 


It obeys the differential equation but violates some of 
the boundary conditions (as is shown in detail in refer- 
ence 5). Hence, Eq. (3) is not a buckling mode, and 
the most that the analysis can yield is that the lowest 
buckling mode can be closely approximated by it. 
The actual deflection function which represents the 
lowest buckling mode can be written as a double power 


series. 


w= D> Dd Omax™y’ (4 


m=0 n=0 


which includes a term representing the simple torsional 
distortion. If it can be shown that this is the dominat- 
ing term, the approximate torsional analysis is justified. 
Instead of the exact deflection function, a close ap- 
proximation to it can be used to demonstrate the pre- 
dominance of the torsional term. This can be con- 
veniently done by the Rayleigh-Ritz Method 


The total potential of a thin rectangular plate without 
any external lateral loads is (see reference 7) 
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For a plate with all edges free, all the boundary con- 
ditions are force and moment conditions. Therefore, 
in the absence of geometric boundary conditions, any 
deflection function is admissible to the minimization 
process of the Rayleigh-Ritz method. If we assume 


a general deflection function 
WwW = 3 >» AnaX” ¥" (6) 


which is a complete set,* the minimization process will 
adjust the coefficients a,,,, so that it represents the actual 
deflection function Eq. (4). If only a finite number of 
terms of Eq. (6) are considered, the minimization will 
make the assumed finite series a close approximation to 
Eq. (4). 

The first 3 terms of Eq. (6), Goo, Giox and doy, are 
rigid body displacements which need not be considered. 


Hence, the first term becomes 
w= ayxy (7) 


which is the same as Eq. (+) and represents the simple 
torsional distortion. 

Substitution of the first k terms of Eq. (6) in Eq. 
5) with the stress distribution of Eq. (2) and minimi- 
zation results in k simultaneous linear equations in the 
coefficients @,,,. These k equations represent the usual 
eigenvalue problem, and the critical temperature param- 
eter /| is given by the lowest eigenvalue. A th 
order determinant is obtained, the lowest root of which 
has to be determined. If successive values for 7},, 
are fairly close to the one obtained for the one-term 
deflection function Eq. (7), the simple torsional distor- 
tion is dominating. 

The stability determinants obtained are in terms of 
constants which depend on the thermal stresses and in 
Hence the effect of aspect 
Increasing aspect ratio 


terms of the aspect ratio.° 
ratio can also be investigated. 
is found to emphasize the dominance of the simple 
torsional distortion, and as the span tends to infinity 
the torsional shape becomes the lowest buckling mode. 

The numerical work is thus aimed at investigating 
the accuracy of the torsional approximation for practical 
cases. An aspect ratio of + has been chosen, since be- 
low that the assumption of infinite wing stress distribu- 
tion becomes unrealistic. The calculations have been 
carried out for solid wings of constant thickness and of 
parabolic cross section, for a number of typical temper- 
ature distributions of the form occurring in accelerated 
flight. The effect of 


variation of aspect ratio between | and 6 within the 


Details are given in reference 5. 


framework of the assumptions of Eq. (2) has also been 


examined. 
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Terms do,y" and da, 9x” of Eq. (6) do not contribute 
to the second integral of Eq. (5), and hence cannot 
represent buckling modes by themselves. They ap 
pear as multiplying constants in the stability deter 
minant having very little effect on the lowest root. 
For symmetrical temperature distributions they do not 
alfect the critical temperature parameter at all (except 
for cases of predominance of a semitorsional buckling 
mode), and for other temperature distributions the ef- 
fect is of the order of 0.02 per cent. Hence these terms 
have been omitted from the analysis. 

The work involved in carrying out the Rayleigh-Ritz 
analysis is reduced considerably if the assumed deflec- 
tion function consists of an orthogonal set, since many 
of the coupling terms will then vanish. In the case of 
the plate it is important that the first and second deriva 
tives of the deflection function also be orthogonal, since 
they appear in the expression for the total potential. 
If such a set can be found, the coupling terms of the 
first integral of Eq. (5) will vanish; unfortunately, the 
stress in the second integral represents a weight func- 
tion which differs from the weight function of the first 
integral, hence the orthogonality in the y direction is 
not preserved and some coupling terms appear. 

In order to represent the deflection function by an 
orthogonal set, it is found convenient to introduce the 


nondimensional coordinates. 


(Ss) 


tr 
I] 
~ 
= 
3 
l| 
2 
= 


For wings of constant thickness the Legendre poly 
nomials very nearly meet the above requirements. 
Since their derivatives are finite series of Legendre 
polynomials of lower degree, the orthogonality of the 
original set is partially preserved for the set of the 
derivatives and only very few coupling terms appear. 
The general deflection function Eq. (6) is therefore ex- 


pressed in Legendre polynomials 
on SE TS aire (9) 
m 1» 1 


where P,,(&) and ?,,(n) are the Legendre polynomials in 
£ and n of degree m and n, respectively. The lowest 
term of the series (9) is the same as Eq. (7) and repre- 
sents the simple torsional distortion. 
In the calculations it is more convenient to compare 
a buckling parameter 
e = akhT, D (10) 


instead of 7}. The successive lowest buckling param- 
eters e, have been computed for three typical tem- 
perature distributions in steel wings of constant thick- 
ness. 

Case 1 (Fig. 3): Instantaneous acceleration to a 
Mach Number of 3 at an altitude of 50,000 ft., the 
boundary layer being assumed to be laminar (see Ap- 


pendix A of reference 6) 


7 = T, — 7,{0.6394 + 0.08940 » — 
0.001074 (1 — n)? (11) 


where 7, = 992.8°R. and 7), = 600.4°R 
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Case 2? (Fig. 3): The symmetric part of Case 1. where Do = Eho*/12(1 — v?) 17 
Case 3 (Fig. 4): A hypothetical temperature distri- 
bution of the type obtained experimentally by rapid 


=> 


Also for the parabolic cross section the successiy 
lowest buckling parameters have been computed for 3 


heating of leading and trailing edges.” ; en ee oes ; j 
typical temperature distributions in steel wings (4 
T = To + T1n' (12) Case 4 (Fig. 6): Instantaneous acceleration to y 
2 ; ; ; ; Mach Number of 3 at 50,000 ft., with a laminar bound. 4 
rhe results of the computations, for an aspect ratio : gt : 4 
; Se ie ary layer (see Appendix B of reference 5) 4 
of 4, are tabulated in Table 1. =e i 
It should be noted that not all & terms of the deflec- ti 7, — 7;[0.9634 (1 — n®) + 1.250(n* — n 
tion function influence the result, since only antisym- 0.4465 (n9 — n>? I8 
metrical terms affect an antisymmetrical buckling : eC 
. 5 where 7\, = 992.S°R. and 7; = 600.4°R. Fic. 2 
mode aid vice versa. However, as the symmetry : ere : eal : eo 
: ; : Case 5 (Fig. 6): The symmetric part of Case 4. 
properties of the lowest buckling mode are not known : Fk 20 ' 
; Case 6 (Fig. +): The same hypothetical temperatur " 
in advance all k terms appear. eee ee ; 
a : , ; : : distribution as Case 3. 
lhe predominance of the simple torsional distortion = : 
: ' rhe results of the computations for an aspect ratio 
is clearly demonstrated since the maximum reduction P ps shia ee 
; ; : pueage ? adi of 4 are tabulated in Table 2. As in Table 1, some of é 
in eis less than | per cent. Case 5 is interesting. The : Sgt 
: ; ; ; the ternis of the deflection function do not influence the 
reduction there is so small that even for a nine-term . ; 
—— Ba a ‘ : result on account of symmetry properties. 
deflection function it is lost in the accuracy of the calcu- : = J : 
lias - Here, in the case of the wing of parabolic cross sec- ¢ 
ations. P : ; , mee : : 
’ ; ; ™ ; one ; ; tion, the simple torsional distortion is not as predomi- is 
Investigation of the effect of variation of aspect ratio ; : . ; ; 
: ‘ . ee ; nant as in the case of the plate of constant thickness. a 
for Cases 2 and 3 indicates very small changes of the ; : f wi 4 
P . a Che shape of the cross section results in reduced stiff 
order of 0.01 per cent even for an aspect ratio of 1. ‘ we 
: : : ‘ j ee ness near the leading and trailing edges, and some bend- 
For wings of parabolic cross section (see Fig. 5) the : Bot d a 
‘ eee ing occurs in addition to the torsional deformation. 
thickness variation is given by : : ; : 
‘ Hence the effect of the higher-order terms is noticeable, é 
h = ho[1 — (y 5)?] (13) and the error involved in neglecting them becomes 
. , ae > slightly larger. Furthermore, a ‘“‘semitorsional’’ mode 
and a particular type of Jacobi’s polynomials, the ultra- ii PEt igs 
‘ : ca : “ (symmetric in x and antisymmetric in y; 1.e., terms of 
spherical polynomials, are a suitable orthogonal set for , aL . 
; i : : Es the type dax°y + de3x*y® + ...) may become the lowest 
the chordwise component of the deflection function. a ae ls : 
é Ag , . buckling mode in certain cases when the compressive ' 
Eq. (6) is therefore written ; 
stresses are concentrated very near the leading and stant 
x we " <a P laver 
; > ¢z ie trailing edges (see Cases + and 5 in Table 2). But as : 
-~ = > »~ Giant m(&)Qn ; (n) (14) . ¥ ad ; ‘ ‘ 
m=1 m=1 the lowest buckling parameter of the “‘semitorsional 
>t ait oa mode is of the same order of magnitude as that of the Vi 
where P,,(&) are Legendre polynomiais in & of degree m . : é ae 
riage : , ; a torsional mode, the error involved in neglecting it is not 
and Q,°: *(m) are ultraspherical polynomials in 7 of I 
is arge. 
order 2 and degree n. The lowest term in series (14) is pis : ; : ; 
rhe analysis was carried out up to terms of the third , 
w= ayP(&) Qi *'(n) dy (15) degree in x and y, and it seems reasonable to exvect on 
ee , = ; ie yhysical grounds that the higher powers, which do not 
which is the same as Eq. (7) except for a multiplying Pht li te | : 3 | 
5 a . . . iy give rise to any new types of shapes, will not affect the 
constant. The buckling parameter is defined for a ez, : ‘ « 
. , . . results considerably. It can, therefore, be concluded ° 
wing of a parabolic cross section as : ie ‘ oe aes ‘ 
that the simple torsional distortion will yield a fair ap- 
e = akh(1\/ Do) (16) proximation to the lowest buckling parameter, the 
TABLE 1 
Lowest Buckling Parameters e for a Plate of Constant Thickness of Aspect Ratio 4 
Lowest Buckling Parameter e, 
Number Case | Case 2 Case 3 
of /) w// of 
Terms below below below 
k Deflection Function Ch e ek e e é 
I an P(E) Pin) 70.022 70.022 22.050 sant 
} anPlE)PilqQ) + ayeP(&) Pon) + ay PE) Pim) 4 69. 580 0.63 70.022 0 22 O50 0 
(to P(&) Poly) 
5 an Py() Pilg) + aeP(é) Ply) + arsP(t) Ply) 4 69572 0.64 70.016 0.01 22.059) 0 
ay P(E) Pi(n) + doxPo(&) Po(n) 
6 a P(E) Pin) + aeP(E) Poly) + aysPi(E) P3(n) + 69.558 0.66 70.016 0.01 22.050 0 
ay PE) P\(n) + dexP2(E) Palm) + ay P(E) Pin 
9g ayP(E) Pin) + aeP(E) Poly) + ayzP(E)P3(n) 4 70.01 0.02 22.050 0 


ay P(E) Pin) + do2P2(E) Palm) 4 
(t93P(£) P3(n) + an P3(€) Pin) + 
AyeP3() Paln) + a33P3() P2(n) 
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TABLE 2 
Lowest Buckling Parameters e for a Wing of Parabolic Cross Section of Aspect Ratio 4 
g ] 
Num- Lowest Buckling Parameter e; . 
ber ——Case 4 Case 5 Case 6 
of w/ % 
Terms below below elow 
k Deflection Function e; é1 e; e e; Pi 
l anP(~é) Qin) 95.605 95.605 46.200 
} anuP(E) Qin) + awP(E) Qo(n) + anP2l(E) Qin) + 94.792 0.85 95.605 0 46.200 0) 
zo P() Q2() 
5 an P(E) Qin) + aywP(E) Qeln) + aisPi(E) Os(n) + 93 .625 2.07 94.055 1.62 15.521 1.47 
ay P(E) Qin) + d22P2() Q2(n) 
6 anP(—) Qin) + awP(E) Qoln) + aisPi(E) O3(n) + 92.755 2.98 93 .920 1.76 45.460 1.60 
ay PE) Qin) + de2P(&) Qe(n) + anPs(E) Qi(n) 
7 anP(~) Qin) + aePi(E) Qaln) + aisPi(E) Q3(n) + 87 .848* 8.11 89.79* 6.08 15.460 1.60 
ay P(t) Qin) + ae2P2(é) Qox(n) + 
d3P(E) Q3(n) + anPs(E) Qi(n) 
g anPy(&) Qin) + awPi(t) Qeln) + aysP1(E) Q:(n) + 89 .79* 6.08 $3.9 $98 
dy P(E) Qin) + dexe2P(£) Qo(n) + (90.4 for tor- 


do3P(E) Qs(n) + as P3(E) Om) + 
432P3(£) Qo(m) + a33P3(€) s(n) 


* “Semitorsional’’ mode predominates. 


error being less than 10 per cent, usually of the order 
of 5 per cent. In practice, in view of the approximate 
nature of the aerothermodynamic data, this error is not 
serious. If greater accuracy is required, a deflection 
function containing many terms must be used. 

Investigation of the effect of variation of aspect ratio 
for Case 6 indicates changes of less than 4.5 per cent 
in the lowest buckling parameter for aspect ratios 
above 2.5. For lower aspect ratios the torsional mode 
is not the lowest buckling mode for wings of parabolic 
cross section. 


DISTRIBUTION OF THERMAL STRESSES IN FINITE WINGS 


In the above buckling analysis the end effect was 
neglected. For wings of small aspect ratio the change 
in the stress distribution due to the free ends may be 
important and has to be investigated. 

The elementary solution for the thermal stress dis- 
tribution, Eq. (2), is exact only for an infinite aspect 
ratio. 
found for wings of constant thickness in the form of an 


For small aspect ratios, an exact solution can be 


infinite series, but the numerical evaluation is difficult 
except in special cases. An example of such a special 
case can be found in reference 9. 

An approximate solution, using the principle of the 
minimum of the complementary energy, was developed 
by Heldenfels and Roberts" for wings of constant thick- 
ness. Another approximate method, developed by 
Horvay'! for the solution of the stress problem of a semi- 
infinite strip, can also be applied here. The two meth- 
ods have been compared in reference 6 and yielded close 
results (see, for instance, Fig. 7 showing the spanwise 
distribution of o, for a plate of constant thickness and 
aspect ratio 3, subjected to the temperature distribu- 
tion of Fig. 3). Horvay’s method of self-equilibrating 
polynomials is more convenient in general, because of 
the complete tabulation of the required functions,'* 
but since in buckling calculations the stresses are re- 
quired in functional form, the first method, which 


sional mode) 


yields the stress distribution as a combination of hyper- 
bolic and circular functions (see Eqs. 1.28—1.30 of 
reference 6), is preferable here. 

Now if the thickness of a solid wing varies, a realistic 
The 


thermal stress problem of a wing of variable thickness 


analysis must take this variation into account. 


differs, therefore, from that of a wing of constant thick- 
ness. However, if the problem is considered in terms 
of stress resultants, rather than stresses, approximate 
solutions can be found for a wing of chordwise varying 
thickness by approaches similar to those of Heldenfels 
and Roberts" and Horvay,'! respectively. A_ stress 
resultant function VY, analogous to the Airy stress func- 
tion, can be defined and the complementary energy be 
expressed in terms of it, and minimized. In reference 6 
approximate solutions are derived by both approaches 
and compared. Another approximate solution, by a 
collocation procedure, is presented in reference 13. 
But again, the method of the minimum of the comple- 
mentary energy yields results in a form most suitable 
Eqs. 1.103-1.105. of 


for buckling calculations (see 


reference 6). 


BUCKLING ANALYSIS FOR WINGS OF SMALL ASPECT 
RATIO 


It has been shown in reference 6, that the thermal 
Stress distribution changes noticeably in the neighbor- 
hood of the free ends x = +a. For wings of small as- 
pect ratio this neighborhood extends over most of the 
wing, and therefore the stress distribution of the finite 
wing has to be used in the buckling analysis. 

The conclusion obtained from the analysis of wings 
subjected to thermal stresses corresponding to an in- 
finite aspect ratio that the simple torsional distortion 
is a good approximation to the lowest buckling mode 
can be expected to hold also for wings of small aspect 


ratio (but above 2.5). 


However, its validity is briefly 
investigated. 


First, it is apparent that the simple torsional distor- 
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3), is also not an exact buckling mode for a 
It can be shown that here, in addition to 


tion, Eq 
short wing. 
the violation of the boundary conditions, also the dif- 
ferential equation is not satisfied. The aim is thus 
again to show the predominance of the torsional distor- 
tion in the general deflection function. 

An analysis along the same lines as before, but using 
the stress distribution for a finite wing (see Eqs. 1.2S 
1.30 of reference 6) instead of Eq. (2), is carried out for 
a wing of constant thickness. Because of the second- 
ary stresses o, and 7,,, the terms a,x” and do,y" of 
Eq. (6 

The successive lowest buckling parameters e,, de- 
fined by Eq. (10), have been computed for the tempera- 


cannot be discarded here. 


ture distribution of Case 1 above (see Fig. 3). The 
corresponding stress distribution for an aspect ratio 3 
is calculated in reference 6, and the results of the buck- 
ling analysis are given in Table 3. 

The calculations indicate a maximum reduction of 
the buckling parameter of the order of 3 per cent, 
slightly larger than in the case of the plate subjected 
to the infinite wing stress distribution, but still small 
enough to verify the predominance of the torsional dis- 
tortion for a short plate of constant thickness. 

It appears, therefore, that the extension of the con- 
clusions about the lowest buckling mode to wings of 
small aspect ratio is valid, though the accuracy of the 
1 


approximation decreases. 


THE END EFFECT IN THE APPROXIMATE TORSIONAL 
ANALYSIS FOR WINGS OF CONSTANT AND VARIABLE 
THICKNESS 


It has been shown that an analysis based on the as- 
sumption that the simple torsional distortion predom- 
inates in the lowest buckling mode, yields a good ap- 
proximation to the critical temperature. This ap- 
proximate analysis is now modified to include short 
wings also. 

If a one-term deflection function kxy, [see Eq. (3) ], 
representing a simple torsional distortion, is substi- 
tuted in the total potential Eq. (5), and the latter is 


minimized with respect to the coefficient k, one obtains 


- D2(1 — v) dxdy + ET, X 


| | . hly?(o,/a@ET\) + 
« be a 


v*(a, a@HT)) + xy(27,,/ak7T)) | dxdy = Q (19) 


In accordance with reference 6, the stresses and stress 


resultants are denoted 


N, = h Cc; = ak Tf "9 ] 
N, = he, = eT, fe" (20) 
N,, = hrr, =. —akl, rg" 


Where / and g are functions of y and x only, respectively, 
whose product defines the stress function or stress re- 
sultant function corresponding to the stress distribu- 
tion of the finite wing. The primes indicate differen- 


tiation of g and f with respect to their assigned inde- 
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pendent variables. For a plate of constant thickness 
Eqs. (20) reduce to Eqs. 1.6 of reference 6, multiplied 
by a constant (a/7) h), and for a plate of chordwise 
varying thickness they reduce to Eqs. 1.93 of the same 
reference, multiplied by a different constant a7}. 
With the assumption that D D(y) and h h(y). 
the substitution of Eqs. (20) in Eq. (19) yields, for a 


unit semi-chord, 


| 
1(1 — v)a | D dy + akT\[A2 C 
* ] 
Fy I] | Ly N, 0 (21 


where 


Fy = } ay 4 
« l 
I] v2 9” dx 29 


Integration by parts of Eq. (22) yields 


l 
- © le l 
« l 


since /(+1) = /’(+1) = O from the boundary condi 


tions at the long edges 


ak Ty fg" = U0 
at y t+) = +] 
~akT\f'g’ = 0 


TABLE 3 
Lowest Buckling Parameters e for a Short 
Plate of Constant Thickness and Aspect Ratio 3 


Lowest Buckling 
Number Parameter ¢; 


of c 


Terms below 
k Deflection Function 
] ay P(E)Pin 89.45 
3 awePAEIPAn) + anP(EIPy) 4 89.43 0 
Ax» PAE)Po(n 
6 aweP A E)P An an P(e)Piln S899 0.49 
ayeP( r U] P. P n 
+ dy PAEVPi(n) + deeP(E)P(n 
7 aoe Po E)P2(y + an P(E)Pi(n S694 2.74 
1eP(E)Ply 43P(£)P3(y 
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Therefore, 


Integration by parts of Eq. (25) yields 


WwW 
= 
7S 
4. 
bo 


since g(+a) = g’(+a) = 0 from the boundary condi- 
tion at the short edges. 


N, = @ilyf'¢ = 0 
at x= +a 
WV os = —akT, fg’ = (0) 
Therefore, He = 3Co (29) 


Similarly, integration by parts of Eqs. (26) and (27) 
yields 
Ll, = A, (30) 
N, == —C (31) 


Substitution of Eqs. (28) to (31) in Eq. (21) gives 
71 
11 — va } Ddy + aET,A2G% =0 (382) 
. l 


since 


Fy H, + Ly Ny = Az Cy — Ao Cy = 0. 


From Eq. (32) it is seen that the effects of the normal 
chordwise stress and the shear stress cancel out. The 
end effect is therefore determined by Cy only, which 
represents the spanwise variation of o;. 

This conclusion holds also for a more general tor- 
sional distortion 


w= G(x) y (335) 


where G(x) is an arbitrary function of the spanwise 
coordinate x. This can be easily shown by repeating 
the above derivation with the more general deflection 
function Eq. (33) instead of Eq. (3). Thus, all the 
conclusions about the end effect arrived at below are 
valid also when the simple torsional analysis approxi- 
mates a case where a “‘semitorsional’’ mode is the lowest 
buckling mode (as, for instance, in Cases 4+ and 5 of 
Table 2). 


If an ‘‘end effect factor’’ is defined by 
/ a 
A = 2a / / gdx = 2a/Co (34) 
/ e a 
the critical temperature parameter becomes 
“1 
Ter = —AZ(1 — ” | D dy aKA, (35) 
I 
now, since the torsional stiffness for thin solid wings is 
+I 
GC = 2(1 — ») / D dy (36) 
« 1 


and, for an infinite span, Eq. (20) becomes 


akT\f" = ho,. = h ET,F(y) (37) 


Aso = 2F (28) 


| gay = 2 | g dx 
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TABLE 4 
The End Effect Factor in the Torsional Analy 


-End Effect Factor \ 


Aspect Plate of Const Parabolic Double 
Ratio Thickness Cross Section Wedgi 
Case 1 Case 3 Case 5 Case 6 Case 7 

] 2 900 2 906 2.393 2.167 2 O82 

2 1.483 1.488 1.401 1. 364 1.349 
3 b.zed 1.280 1.233 Lae 1.209 

} 1.194 1.196 1.167 1.154 1 149 

5 1.150 1.15) 1.128 1.120 1.115 
6 1. 122 1.123 1.105 1.098 1.094 
7 1.108 1.108 1.089 1.083 1.080 
8 1.089 O89 1.076 1.072 1.069 


Eq. (35) can also be written, for comparison with Eq 
(1), as 


*! 
Tice = DX —GC/ak | hy? F(y) ay | OS 
1 


Eq. (38) is simply the result of the torsional analysis 
multiplied by the end effect factor. \ can thus be re- 
defined as the ratio between 7),, for a plate of finite 
aspect ratio to that for a corresponding infinite plate. 

Since g is given for a wing of constant or variable 
thickness by Eq. 1.22 of reference 6, \ can be calculated 
from Eq. (34), once the approximate solution to the 
thermal stress problem is known. From Eq. 1.22 of 
reference 6 write 


Co = | gdx = 2a + [2 (p* + q*)|[(Capb + Cig) X 
cosh pa sin ga + (Cip — Cig) sinh pa cos ga 39 
By substitution for C; and C; from Eqs. 1.25 and 1.27 
of reference 6, Eq. (39) can be written as 
Ce = Be — [4pg (p- + q°) x 
cosh? pa sin? ga + sinh? pa cos? ga | 
1() 


p cos ga sin ga + g sinh pa cosh pa 
Now, for pa > 5, cosh pa = sinh pa 
and cosh pa sinh pa > cos ga sin ga 


Then Eq. (40) reduces to 


Cy = 2a — [4p /(p? + g*)] t] 
For a> om, Co 2a 
and, from Eq. (34), A= I 


Values of \ have been computed for various aspect 
ratios for the temperature distributions of Cases | and 
3 of the constant thickness plate, for Cases 5 and 6 oi 
the wing of parabolic cross section, and for Case 7, a 
double wedge section subjected to the same tempera- 
ture distribution as Cases 3 and 6 (see Fig. 3). The 
results are tabulated in Table 4 and plotted in Fig. § 
The values of \ for Cases 1 and 3 of the constant thick- 
ness plate are very close; similarly for Cases 5 and 6 
of the wing of parabolic cross section. For this reason 
only Cases 3, 6, and 7 are drawn in Fig. 8. It appears 
that the magnitude of the end effect is not changed con- 
siderably by the temperature distribution, since this 


(Continued on page 590) 
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Flight Mechanics and Variational Problems of 


a Linear Type’ 


ANGELO MIELE* 


Douglas Aircraft Company, Inc. 


SUMMAR\ 


4 general analysis is presented of variational problems of 
linear type. The characteristic of these problems is that the 
functional forms to be extremized, the differential equations of 
constraint, and/or the eventual isoperimetric conditions only in 
yolve the first power of the derivatives of the unknown functions 
For the above problems the Legendre-Clebsch condition fails to 


vield any information on the minimal or maximal nature of the 


moreover, the Weierstrassian function is zero along 


extremal arc; 
By the use of a new 


, major portion of the Eulerian solution 


technique, based on Green’s theorem, the afore-mentioned 


difficulties can be overcome. Rigorous sufficiency proofs are 
developed for linear problems of both the simplest type and the 
isoperimetric type. Several physical problems, associated with 
the mechanics of flight of either a constant mass or a variable mass 
vehicle, are treated within the general frame of problems of 


linear ty pe 
(1) INTRODUCTION 


S VERAL NEW PROBLEMS of applied mathematics 
have arisen in the analysis of trajectories of high- 
speed aircraft and missiles which cannot be handled by 
the use of conventional methods of performance anal- 
ysis. These problems, which fall under the domain 
of the calculus of variations, are generally concerned 
with the determination of the maneuver of the em- 
pennage control and or the engine control which opti- 
mize a given performance. 

A category of variational questions of rather fre- 
quent occurrence in aeronautical engineering is repre- 
sented by problems of the so-called linear type. The 
dominant characteristic of such problems lies in the 
fact that the functional forms to be extremized, the 
differential equations of constraint, and or the eventual 
isoperimetric conditions only involve the first power of 
the derivatives of the unknown functions. This ap- 
parently simplifying circumstance is actually the source 
of important analytical difficulties. In connection with 
the study of the minimal or maximal properties of 
the solutions, the Legendre-Clebsch condition fails to 
vield any information on the nature of the extremal arc. 
An analogous difficulty is found, along a major portion 
of the extremal arc, when applying the more stringent 
Weierstrass condition. 

$y the use of a new and hitherto unexploited tech- 
nique, based on Green's theorem, this difficulty can be 

Received February 13, 1958. Revised and received April 
18, 1958 

rhis article is a condensed version of the investigation de 
scribed in reference 11 
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overcome. Rigorous sufficiency proofs can be devel 
oped for linear problems of both the simplest type and 


the isoperimetric type, as is shown below. 
(2) LINEAR PROBLEMS OF THE SIMPLEST TYP! 


Consider a functional form of the typet 


Ox 
[ = | [P(x, vy) + v(x, vy) dx | 
where ® and W are given functions of the arguments a 
and y and where y’ = dy dx. Assume that, among all 
arcs v(x) consistent with certain prescribed end condi 
tions, the special are which extremizes—1i.e., maximizes 


the integral (1) is to be found. As is 


or minimizes 
known from the theory of the first variation, the ex 
tremal are associated with the problem under consider 


ation is to be a solution of the Euler's differential equa 
tion 


(d dx) (OF, Oy’ (OF oy 0 2 


where the so-called fundamental function F is defined as 
F = P ‘i Vy’ BS | 


The development of Eqs. (2) and (3) leads to the fol 


lowing result: 


(OV Ox) — (0% Oy) = 0 | 

The above equation points out that there is, asso 
ciated with the functional form (1), a degenerate vari 
In fact, while the Euler's differential 


ational problem. 
differential 


equation is, ordinarily, a 
equation (as such, its general integral depends on two 
(4), on the contrary, is an 


second-order 


integration constants), Eq. 
equation in finite terms. 

It is, therefore, concluded that if the coordinates of 
the end points | and 2 are consistent with Eq. (4), then 
Eq. (4) represents the extremal are associated with the 
variational problem under consideration. If, on the 
contrary, the end points 1 and 2 do not lie on the curve 
defined by Eq. (4), then the variational problem has no 
solution. 

It is to be noted that, even when the boundary con 
ditions of the problem are consistent with Eq. (4), a 
further difficulty arises. For a linear integrand the 
Legendre function 0°F Oy” is zero at all points of the 
the Weierstrassian / func- 
It is, there 


extremal arc. Moreover, 
tion! is also EH = OQ along the extremal arc. 
t The subscripts 1 and 2 are employed here to denote, respec 


tively, initial point and final point 
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Fic. 1. Parametric representation of a variable subject to 


inequality. 


fore, inferred that both the Legendre condition and the 
Weierstrass condition yield no information on the 
minimal or maximal nature of the extremal path. 


(3) LINEAR PROBLEMS INVOLVING A NONHOLONOMIC 
CONSTRAINT 


With regard to engineering applications, a more im- 
portant type of linear problem occurs in the case where 
the are y(v) extremizing the integral (1) is to be con- 
sistent with a constraint of the form 


where z is a new variable such that* 
Ss 5S & (6) 


To account for the relation (6), the idea of para- 
metric representation of the variable subject to in- 
equality is introduced. The variable z is represented 
(see Fig. 1) as a function 2(p) of a parameter p having 
the following properties: (a) for —®e < p < py, 
z= %, dz/dp = 0, d’z/dp? = 0; (b) for pf < p < 
+ ©, = 20, dz/dp = 0, d’z/dp? = 0; and (c) for p) < 
pb < ps, a one-to-one correspondence is assumed between 
sand p—i.e.,dz dp ¥ O everywhere. 

With the above artifice, the independent parameter 
p is allowed to vary between —o and +o. The 
variable z becomes a dependent quantity, variable be- 
tween 2, and 2. according to the scheme of Fig. 1. 

Notice that the condition dz dp = 0 represents either 
a line zs = 2, ora line zs = zg; on the other hand, the con- 
dition dz dp ¥ 0 is representative of any other value 
of z intermediate between 2, and z». Notice also that 
p is only a parameter and that there is no necessity of 
attributing to it any geometrical or physical meaning. 

The main effect of the parametric representation is 
to reduce a variational problem involving inequalities 
to the same mathematical model useful in solving the 
case where all constraints involve equalities. Thus, 
for the present problem, the two unknown functions 
y(x) and p(x) must be determined so as to extremize 


* The subscripts 1 and 2 are employed here to denote, respec- 


tively, lower bound and upper bound for the s variable. 
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the integral (1) under the new constraint 
y’ + flx, ¥, 2(p)] = 0 7 
The dependent variables being two in number, ty 
Euler equations are needed to define the extremal ar 
(d dx) (OF Oy’) — (OF oy) = 0 g 
(d dx) (OF Op’) — (OF Oop) = 0 (9 
F being the so-called fundamental function 
F= 6+ Wy’ + Xy' +f) (10 
and A(x) a variable Lagrange multiplier. Eqs. (8) to 
(10) lead to the following results: 
(d dx) (W + Xr) — (0@ Oy) — 
(OW Oy)y’ — A(Of Ov) = 0 (11 
(Of 0s) A(dz dp) = 0 (12 
The Euler equation (12) associated with the /p variable 
is interesting, because it shows that the extremal arc is 
discontinuous, being generally composed, for Of Oz ¥ 0 
of subares 
A = 0 (13 
dz/dp = 0 (14 


and subares 


For the subare \ = 0, the relationship between » and 
y is obtained from Eq. (11) 


(OW Ox) — (0b dy) = 0 (15 
For the subarcs dz/dp = 0 the z variable has a con- 
stant value—i.e., either 2 = 2, or 3 = 2» according to 


the parametric representation devised in Fig. 1. The 
relationship between x and y is obtained by integrating 
Eq. (5) subject to either the constraint s = 3, or s = 2. 

Clearly, the present problem differs from the one of 
Section (2), insofar as subares s = const. must be con- 
sidered in the composition of the extremal are in addi- 
tion to the special subare represented by Eq. (15). The 
role of these subares 2 = const. is fundamental in the 
solution of the boundary value problem: in fact, it is 
now entirely possible to satisfy a set of prescribed end 
conditions in the (xv, y) plane—a circumstance which, 
on the contrary, is generally impossible for the problem 
formulated in Section 2. 

The analysis of the Legendre-Clebsch condition and 
of the Weierstrass condition (omitted for the sake of 
brevity) leads to the following results:'! (a) the 
Legendre-Clebsch condition yields no information on the 
minimal or maximal character of the solution for the 
three subarcs A = 0, 2 = 2, and gz = 22; and (b) the 
Weierstrassian function is / = 0 at points of the sub- 
arch = 0. The analogy between the present conclusions 
and those of Section (2) should be emphasized. 


(4) LINEAR PROBLEMS INVOLVING BOTH A 
NONHOLONOMIC CONSTRAINT AND AN ISOPERIMETRIC 
CONSTRAINT 


As a further step, the problem of extremizing the 
functional (1) is considered in connection with a 
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FLIGHT MECHANICS AND 
nonholonomic constraint obeying Eq. (7) plus an iso 
perimetric constraint of linear type: 


xX 
| } [P(x, vy) + W(x, yy’ |dx = given (16) 
J x 
For the problem represented by Eqs. (1), (7), and (16) 
the two unknown functions y(x) and p(x) must be con- 
sistent with the Euler equations (S) and (9) where the 


F function is now defined as follows: 
F= 6+ Wy +Xy' +f) +A + Wy’) (17) 
\(x) being a variable Lagrange multiplier and \ a con- 


stant Lagrange multiplier. The development of Eqs. 


M(Of Oy) — A[(OH/Oy) + (OW Oy)y’] = O (1S) 


S), (9), and (17) leads to the following result: 
d dx) (VW + XW +A] — (08 dy) — (OW dy)y’ — 


(Of Oz) A(dz dp) = O (19) 


Once more the extremal are is discontinuous, being 
composed, for Of 02 ¥ 0, of subares dz dp = 0 (s = 2, 0r 


¢ = 2) and subarcs \ = O—i.e., subares such that 
(OW, Ox) — (by dy) = 0 (20) 
where Vv, = V+ iv (21) 
?, = b+ (22) 


Notice that Eq. (20) is representative of a family of 
curves in the (v, y) plane, one for each value of the con- 
stant multiplier X. The particular value of \ asso- 
ciated with the problem under consideration must be 
computed on the basis of the given isoperimetric condi- 
tion and of the prescribed boundary conditions. 

The study of the Legendre-Clebsch condition and of 
the Weierstrass condition is omitted here, as it leads"! 
to conclusions similar to those of Sections (2) and (3). 


5) TypIcAL PROBLEMS INVOLVING A NONHOLONOMIC 
CONSTRAINT 


In the present section several physical problems are 
introduced; their mathematical model falls within the 
frame of the problem considered in Section (3). 


(5.1) Burning Program for the Maximum Range of a 
Rocket-Powered Aircraft in Level Flight 


A rocket-powered aircraft moving along a rectilinear 
horizontal path is considered in connection with the 
a) the thrust is tangent to the 


following hypotheses : 
flight path, (b) the equivalent exit velocity of the 
rocket engine (I”,) is a constant, and (c) the engine is 
capable of delivering all mass flows (8) bounded be- 
tween a lower value (8 = 0) and an upper value 
] Bb, 

In flight at constant altitude, the drag D can be con- 
ceived as a function D(I’, £) of velocity (V’) and lift 
L) only. After accounting for the equation of motion 
on the normal to the flight path (1 — m g = O, where 
m is the mass and g the acceleration of gravity), the 
drag D is reduced to a function D(V, m) of velocity 


and instantaneous mass of the vehicle. 
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The horizontal distance (XY) flown by the aircraft is 


given by* !* 


; ai ml" VV 
A = | ~ -—— —m' idV (23 
Ji Dim, V) Dim, V) 


where m’ = dm dV. The maximum range problem 
consists of extremizing the integral (23), subject to the 


constraint 
m' + j}meB [VB — Dim, V)]}{ = 0 (24 


Baz (25 


lA 


where 0< 8B 


It is to be noted that the problem represented by the 
relations (23), (24), and (25) is mathematically identical 


with the problem represented by Eqs. (1), (5), and (6). 
This fact is evident if the following substitutions are 
made: 
x=V, y=m, 2=8, 3%=0, 2 B, 

® = —(mV D), VW = —-(I,V D), 

f= mBp'(V.B— D 
The solution, therefore, involves (a) subares 8 = 0, (b 
subares 8 = Byer, and (c) subares flown with variable 
thrust, along which the mass-velocity relationship is 


represented by 
(0/OV’) [VV D] — (0/0m) [mV D] = 0 26 


For the particular case of a drag polar of the form 


D = AV? + B(m?/V 2% 

(where A and B are appropriate constants), Eq. (26 
reduces to!” 

m= VA/BV?V(V4+1),(V+3V) (28 


With regard to the boundary value problem, the 
subares 8 = O and @ = 8B,,,; must be combined with 
those obeying Eq. (26), so as to satisfy a set of pre 
scribed end conditions. In this connection, the further 


inequality 
dm <0 (29) 


Its effect is to limit the number 
This 


must be accounted for. 
of physically possible combinations of subarcs. 
circumstance is described in Section (7). 
(5.2) Burning Program for the Maximum Endurance of 

a Rocket-Powered Aircraft in Level Flight 

Under the same hypotheses of Section (5.1), the 
problem of maximum endurance consists in extremizing 
the integral 

V m V, ; a 
t= _ -_—— —~m idV (30 
Jv; Dim, V) Dim, V) 

subject to the constraints (24) and (25). The solution 
involves subares 8 = 0, subares 8 = 8,,,;, and variable- 
thrust subarces satisfying the relationship 


(0 OV) (V..D) — (00m) (m D) = 0 (31) 


For the particular case of a drag polar obeying Eq. (27 
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the variable-thrust subarc reduces to 
m= VA BY? (32) 


(5.3) Burning Program for the Maximum Increase in 

Altitude of a Rocket in Vertical Flight (Zero Drag) 

A rocket-powered vehicle moving along a rectilinear 
vertical path is considered in connection with flight 
in a vacuum. Hypotheses (a), (b), and (c) of Section 
(5.1) are maintained. Assuming that the initial alti- 
tude is zero, the problem of maximum altitude increase 


consists in extremizing the integral‘ 


»\ 
h= } [—(l g) — (VV mg)m’ |dV (3:3) 


| 


subject to the constraint 


m’ + [Bm (IB — mg)] = 0 (34 
where OS BS Bou (35) 
The solution involves only subarecs 8 = O and subarcs 
8 = Bar Variable-thrust subares do not participate 


in the composition of the extremal are insofar as 


(0 Ol’) [V.V mg] — (0 om) [lV g] = 


V. mg #0 (36) 


(5.4) Burning Program for the Maximum Increase in 
Altitude of a Rocket in Vertical Flight (Quadratic Drag 
Law) 

The problem of maximum altitude increase is here 
investigated once more. All the hypotheses of Section 
(5.3) are maintained, the only exception being that a 
quadratic drag law (D = AI”, with A = const.) is now 
assumed. 

The analytical question consists in extremizing the 


integral* 
h= | }—[mV (mg + KIV2)] - 
Jv 
[VV (mg + KV2)]m’} dV (37) 


subject to the constraint 


ae 


m’ + [Bm (V8 — mg — KV?*)] = 0 


where <2 <5 8... (39) 


The solution involves subares 8 = 0, subares B = Bryor. 
and subarcs flown with variable mass flow. Along the 


latter, the mass-velocity relationship is represented by 


(0 OV) [V.V/ (mg + KV?)] - 
(00m) [mV (mg + KV?)] = 0 (40) 


The development of Eq. (40) leads to 
m = (KV?2/g) [1 + (V/V,)] (41) 


(5.5) Brachistocronic Climbing Technique for a Turbojet 
Powered Aircraft 
The climbing flight of a turbojet-powered aircraft 
is now considered in connection with the following hy- 


potheses: (a) the weight (IV) is ideally a constant; 


(b) the equation of motion on the normal to the flight 
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path is approximated as L — II’ = 0, where L is thy 
lift; (c) the thrust (7°) is tangent to the flight path; (q 
the power plant is operated in such a way that thrust 
(7°) and fuel consumed per unit time 
functions 7( 1", h 


(g) are known 
and g(J’, 4) of altitude (4) and y 
locity ( 1’). 

For flight in a vertical plane the drag D can be con 
ceived as a function D(h, \V’, L) of altitude, velocity 
and lift. 
drag D is reduced to a function D(h, V 


After considering hypotheses (a) and (b), thy 
) of altitude and 
velocity only. As a consequence, the time (f) neces 
sary to transfer the aircraft from an initial conditior 
| toa final condition 2 is given by® 6 


(1 — D)| +4 


QV) (W(r-—Dijh'f dV (yp 
where h’ = dh dV. 

The brachistocronic problem consists in finding the 
altitude-velocity relationship which minimizes the in- 
tegral (42) subject to the constraint 


h'+ (Vg) [WV sin 6 (W sing + D — T QO (43 


where @ is the path inclination with respect to a hori- 
zontal plane and 


—l<sinéd < +1 (44 


‘his problem can be reduced to the mathematical mode 
treated in Section (3), by using the substitutions 


x=V, y=h, z=smn0, 3 = -1 2 =4 


@= W/giT —D), vV=Wd(T-D 


f= (I g)[WVsiné (W sing + D— 7T)] 


The solution, consequently, involves (a) subares sin 
6 = —1 (vertical dive), (b) subares sin @ = +1 (vertical 
zoom), and (c) subares flown with variable path incli- 
nation, along which the altitude-velocity relationship 1s 
represented by* 


(00V) [(W VT — D)| - 


(00h) [IV g(T — D)] = 0 (45 


With regard to the boundary value problem, the 
subares sin 6 = —1 and sin @ = +1 must be combined 
with those cbeying Eq. (45) so as to satisfy a set of pre- 
scribed end conditions. In this connection, the case 


where the further inequality® 
T—-D>0 16 


must be accounted for is of engineering importance. 
Its effects are to limit the number of physically possible 
combinations of subarcs. 


* The analysis of the climbing technique of minimum fuel 
consumption® © Jeads to a mathematical model similar to the 


one of Section (5.5). An analogous remark refers to the so-called 


~ 


steepest climb, if the simplifying assumption sin 6 = tan @ is 


accepted. 
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(5.6) A Degenerate Case: Burning Program Extremizing 
the Time of Flight for a Rocket in Vertical Flight (Zero 
Drag) 
With the same hypotheses of Section (5.3) the prob- 
lem of extremizing 


>] 
f(1/e) =~ (V, 
1 


is considered, subject to the constraints (34) and (35 


mg)m'|dV 17 


The solution formally involves subares 8 = 0, sub 
ies 8 = Bar, and subarcs such that 

oO Ol) [I. mg] — (0 Om) [1 g] =O 1S 

Notice that Eq. (48) degenerates into an identity;* this 


isa clear indication of the fact that an infinite number 
f extremals exist,’ all leading from the initial point 

m,) to the final point (I, m2) in equal times. 
[he latter circumstance is logical, insofar as Eq. (47 


can be integrated in a closed form yielding 
ls + V, log (m, me) ] (49 


PROBLEMS INVOLVING BOTH A 
AN ISOPERIMETRIC 


6) TYPICAL 
CONSTRAINT AND 


CONSTRAINT 


NONHOLONOMIC 


In the present section several physical problems are 
considered: they all fall under the mathematical scheme 


of the problems formulated in Section (4). 


(6.1) Burning Program Maximizing the Range of a 
Rocket-Powered Aircraft for Given Time (Level Flight) 
(5.1) 


5.2), the problem of maximizing the integral (23) is 


Under the same hypotheses of Sections and 
considered, in connection with the nonholonomic con- 
straint represented by the relations (2+) and (25) and 
with a further isoperimetric constraint, represented by 
Eq. (30), where ¢ is a given quantity. 

This problem is mathematically identical with the 
one represented by Eqs. (1), (7), and (16), as is evident 
from the substitutions 


fom § y=m, z:=8, 3=0, 2=8 
b=-mV D, V= -V.V D, f= mp (V8 —-— D) 
@’= —-m D, w= —-V,D 

The solution arc, therefore, involves subares 8 = 0, 
subares 8 = By.;, and variable-thrust subarcs. The 
latter are ruled by the follwing relationship: 
001) (1, D) (V+ XJ - 

(0 Om) [(m D) (V+ X)] = 0 (50) 


Notice that Eq. (50) represents a one-parameter family 
of curves in the (m, I’) plane; each member of the 
lamily corresponds to a different value for the constant 
multiplier X. Once the isoperimetric condition and the 
boundary conditions are prescribed, the particular 
value of \ associated with a given problem can be de- 


termined. For the particular case of a drag polar 


* Therefore, any arbitrary function m(V) can be regarded as 
1 solution to Eq. (48). 


\ 
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obeying Eq. (27), Eq. (50) reduces to 
mv+ 0 


wink Vir sar PV +2) _ 


(6.2) Brachistocronic Burning Program of a Rocket in 
Vertical Flight for Given Altitude Increase (Quadratic 
Drag Law) 

Under the same hypotheses of Section (5.4), the prob 


lem of minimizing the integral' 


whe [ m | 
7 I [ me + AT me + KV? tie - 


The’ nonholonomic constraimt 
sented by the relations (3S 


together with the isoperimetric constraint represe! ted 


is considered. repre 


and (39) is accounted for 


by Eq. (37), where /i is now a prescribed quantity 
Q, subares 6 = 3 
the latter 


The solution involves subares 8 


and variable-thrust subares. Concerning 


the mass-velocity relationship is represented by 


(0 O17) | [V. (mg + AV*)] CU + AV); 
(0 Om) | [m (mg + AV2)] 1 + XI 33 

i.e., by 
m= (KV gl.) (VV + Vi) + (UV +210) (O4 
(6.3) Brachistocronic Climbing Maneuver of a Turbojet 

Powered Aircraft for Given Horizontal Distance 

Under the simplifying assumption sin 6 => tan @, and 
conserving all other hypotheses of Section the 
horizontal distance (Y) traveled by a turbojet powered 


aircraft while climbing is given by 


pep ow Ww 1, 
nlLeT-D'T-D | 


The problem of minimizing the integral (42) 1s now con 
sidered, subject to the constraints (43) and (44) and to 
the further isoperimetric constraint represented by 
Eq. (55), where X is a prescribed quantity.’ 
The solution are includes subares sin 6 sub 


+1, and subares along which the eltitude 


ares sin 6 = 
velocity relationship is expressed byt 
- D)\j 

1 — D I (} ob 


WV 
(0 Oh) | (IW g) [U1 


(0 O] (ld + AV) (7 
+ XV 


7) APPLICATION OF GREEN'S THEOREM TO THI 


AND (5 


PROBLEMS OF SECTIONS (3 


In the previous sections a general treatment has been 
developed, valid for linear problems. Several physical 
All of them involve 
the an- 


problems have been formulated. 
discontinuous solutions. For each of them 
alytical relations governing the different subares form- 
ing one extremal arc have been indicated. For all of 
them the Legendre-Clebsch condition fails everywhere 

+ The problem of the climbing technique of minimum fuel 
consumption for given horizontal distance’ leads to a similar 
mathematical model and, therefore, to a solution analogous to 


the one described here 
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and the Weierstrass condition fails along a major por- 
tion of the extremal arc. 

In the present section, Green's theorem is employed 
to overcome the difficulties associated with the indirect 
methods of the calculus of variations in the analysis of 
linear problems. A method of attack is outlined by 
means of particular examples having physical interest. 
The technique here introduced is general. In fact, with 
minor variants, it applies to all questions formulated 


in Sections (3) and (5), and to other problems, also.” ! 


(7.1) Maximum Increase in Altitude for a Rocket in 
Vertical Flight (Zero Drag) 


) 


The problem of Section (5.3) is now considered. 


After introducing the definitions 


@= —(I' g) WwW = —(I,I mg) (57) 
Eq. (33) is rewritten as 
al 
h= | (Pb + Wm’ jd (5S) 
JV 


(7.1.1) The w Function—The following function of 


the variables |” and m is defined: 
w = (OW Ol’) — (0% Om) (59) 
After accounting for Eqs. (57) and (59), one obtains 
w = —(V,./mg) (60) 


The w function, therefore, has the property of being 
negative at all points of the semiplane m > 0. 

(7.1.2) Class of Admissible Displacements in_ the 
Mass-Velocity Plane—Consider a rocket having the 
state of flight (I’, m) at time instant ¢ and let P be its 
representative point in the mass-velocity plane (see Fig. 
2). After a time interval dt = t’ — ¢, the new mass 
of the rocket is m’ = m + dm, the new velocity is 1” = 
1° + dl’, and the new vector position of the rocket in 
the (V, m) planeis OP’ = OP + PP’. 

As is evident from Eq. (34) and from the inequality 
dm < 0, there is one elementary displacement PP’ for 

Since 0 < B < Bar, one 
admissible displacements 


each value of the mass flow @. 
infers that the class of 
PP’ is bound by two limiting conditions 
ment PR(3 = 0), corresponding to coasting flight; and 
a displacement P7’ (8 = Bn.,), corresponding -to flight 
with maximum engine output. 

While a displacement like PS is physically possible, a 
displacement like PQ is physically impossible, because 
the mass of the rocket cannot increase; analogously, a 
displacement like PZ is also physically impossible, be- 
cause the constancy of the mass implies zero thrust and, 


a displace- 


therefore, a negative acceleration. In conclusion, the 
two lines 8 = Oand 8 = Byq, split the (I, m) plane into 


a region (//), 


two regions associated with point P- 
whose points are accessible to the rocket in vertical 
flight; and a region (A), which is forbidden to the 
rocket, because of the physics of the motion. 

(7.1.3) Analytical Proof 
must be transferred from the initial state of flight 
(1), m) to the final state of flight (I, mz»). 


Consider a rocket which 


As a first 
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step, the two limiting lines 8 = O and 8 = 8,,q, starting 
from | are calculated. Their respective equations ay 


” 


(see Fig. 3) 
m = const. (()] 
Io + I, log m — (gm Bnar) = const. (j2 


As a second step the two other limiting lines 3 = 
and 6 = By, ., arriving at 2 are computed (see Fig. 3 
The broken line 13241 of Fig. 3 is drawn: it surrounds 
a closed region to which all admissible paths belong. 

Now focus attention on the Eulerian are 132 whic! 
is formed by one subare 13 flown with maximum en 
gine output followed by one subare 32 flown by coast 
ing. The arbitrary, though admissible, compariso 
path 152 is considered and the difference in final alt 
tude between the paths 132 and 152 calculated 


Ah = hye — Is. = (b + Wm’) dV — 
7! 


/ (db + Wm')dV = g (ddl + Wdm) (638 
¢ 152 13251 


By means of Green's theorem, the above cyclic integral 
can be transformed into a surface integral associated 
with the area a (see Fig. 3) bounded by the clockwisy 


circuit 13251: 


Ov OP : 
a= — ff _— dVdm = 
JJa \Ol om 
— | wlldm (64 


/JdJ/a 


Since the w function, defined by Eq. (59), is negative 
everywhere, it follows that 


Ah > 0 (65 


In view of the arbitrariness of the comparison path 152 
it is concluded that the Eulerian are 132 maximizes 
the altitude rise, under the limitation 0 < 8B < B,,,,. 

By the use of the same technique, it can be proved 
that the other Eulerian path 142 (see Fig. 3), formed 
by a coasting subare 14 followed by a subare 42 flown 
with maximum thrust, has the property of minimising 
the altitude rise /. 


(7.2) Maximum Range of a Rocket-Powered Aircraft in 
Level Flight 
The problem of Section (5.1) is now considered 
After introducing the definitions 
@= —(mV/D), ¥ = —(V.V/D) (66 
Eq. (23) is rewritten as 
al” 
X = | (b + Wm')dV 67 
a if 
(7.2.1) The w Function—The following function ol 
the variables I” and m is defined: 


w = (OW Ol’) — (Of Om) (GS 
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I. {1(0OD OV) — DJ + VID — m(OoD om)] 
Ww Dd? 
(69) 


For a drag polar having the form represented by E 


(27), the w function becomes 


AV?(V + V,) — Bim? V2) (V + 3Y, 
a mM 


IAI? + Bim? | 


(70) 


and, therefore, embodies the following properties (see 
Fig. +): (a) w = O, at points of the variable-thrust 
subare, defined by Eq. (28); (b) w < 0, in the region A 
bounded by the m axis and by the subare governed by 
Eq. (28); and (c) w > O, in the region B bounded by 
the |” axis and by the subare governed by Eq. (2S). 

All of the subsequent line of reasoning refers to the 
simplified case represented by Eq. (27)—1.e., the case 
of a drag polar with constant coefficients. The analy- 
sis of a drag polar with coefficients depending on the 
Mach Number has been carried out by the writer in 
reference 2. 

(7.2.2) Class of Admissible Displacements in_ the 
Mass- Velocity Plane—The class of admissible displace 
ments for the rocket-powered aircraft in the mass- 
velocity plane is entirely defined by Eq. (24) and by the 
inequalities (25) and (29) and is determined by con- 


siderations analogous to those of Section (7.1.2) 


In particular, it is of interest to trace the two limiting 
lines 3 Q and g Bmar Starting at | and the other 
two limiting lines @ = O and 6 = £@,.,; arriving at 2; 
the associated broken line 1532641 surrounds a closed 
region to which all admissible paths belong (see Fig. 5 

(7.2.3) Combination of Subarcs—Four main types of 
boundary conditions may exist, depending upon the 
relative position of points | and 2 with respect to the 
curve w = 0. Four cases of flight must be, accordingly, 
analyzed. Generally speaking, the optimum path 
1562 includes three subares: 15, 56, and 62. The 
subare 56 is always flown with variable thrust accord 
ing to Eq. (2S); the subare 15 is flown with maximum 
engine output if point | belongs to the region A or by 
coasting if point | belongs to the region B; the sub- 
are 62 is flown with maximum engine output if point 2 
belongs to the region B or by coasting if point 2 be- 
longs to the region A. 

(7.2.4) Analytical Proof—An analytical proof is here 
offered for one of the four cases of flight listed in Table 
l~-namely, the case of flight designated as Case I. 
For all the other cases, the demonstration is analogous. 

The optimum path 1562 is compared with an arbi- 
trary, though admissible, control path 172 and _ the 


following difference is formed: 
AX = X 1562 = X72 _ 


(bd + Wdm) + (ddl + Wdm) (71) 





1571 7627 


By the use of Green's theorem, the above two cyclic 


integrals are transformed into surface integrals 
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TABLE | where 
Region Region 
where where Extremal Ar 
Case of — Point 1 Point 2 Subare Subare Subar 
Flight belongs belongs 15 56 62 
] \ B B8B=8 wo = 0 = 2 
II \ \ 8B=8 w = 0 ( ‘ 
| a - Sa 
II] B \ B=0 w = 0 the vi 
I\ B B B= 0 wo = 0 
fig iv simp! 
Ax = = wdVdm + | w dVdm 7? 
we SA ee SB 
° . e ° ° . ° A 
of which one is associated with the area S, internal t 
the clockwise contour 1571, the other one with the are 
S, internal te the counterclockwise contour 7627. I; Notic 
Note 
view of the fact that w < 0 at points of the area S, and ' 
. fund 
w > Oat points of the area Sz, it follows that 
AX > 0 73) | are, d 
boun¢ 
Since the comparison path 172 is arbitrary, the in- | py Ec 
equality (73) proves that the extremal are 1562 ma bet we 
mises the range for given end velocity and given pro- On 
pellant mass. the s 
reg1oi 
(S) APPLICATION OF GREEN'S THEOREM TO THI indice: 
PROBLEMS OF SECTIONS (+) AND (6 R.1 
; ae Ply, . mut ti 
By means of minor modifications of the consider- i 
. ‘ ‘ : - ‘ missi! 
ations developed in Section (7), Green's theorem cai ; 
. a . . ill¢ 
also be employed to establish sufficiency proofs in 
; a ' ; : ; Cons 
connection with problems involving constraints of the 
: : ; . must 
isoperimetric type. As an example, the problem of 
‘ : > : : : : j nes 
Section (6.2) is investigated. It should be noticed ; 
; . 1" 2, the 
that the present theory is general and applies to all ; 
c CO : ‘ 1 , i clos 
problems of Section (6), as well as to other problems. Rig 
10, 
(8.1) Brachistocronic Burning Program of a Rocket in Ph 
Vertical Flight for Given Altitude Rise (Quadratic Drag tions 
Law ‘ 
) that 
The problem of Section (6.2) is now investigated subar 
The following definitions are introduced: subar 
= being 
© = —(|m (me + AI), W = ee 
—|l. (mg + AIT*)] (74 re p 
- ' 3 , ees : satis! 
@ = —([ml (me + AI*)], VW = 
—([I.V (mg + AV nso 
i = ticula 
and the integrals (52) and (37) are rewritten as condi 
- III oi 
; = | (®p = Wm’ \dV 76 Col 
JI paris 
+I 
h= } (b + wm')dV (77) § 4 
J} 
After introducing the further definition 
r=t+0h (78 
By m 
one obtains are tr 
(79 


. 
r= [ (dy + Wam')dV 
ea Vi 





nal t 


, and 


pro- 


t in 
Drag 


i) 


76 
‘ 


if 


FLIGHT MECHANICS AND 
where 
dy > + XP (SO 
Vv, = V+ dw (S] 
S1 The wx Function—The following function of 
the variables I’ and m is introduced 
a (OWx O17) — (Os Om (S2 
Simple m inipulations yield the following results 
AV +21.) + AAT7(0 + iV) — Ame I, 2 
s. 


(mg + AI)? 


Notice that, for a given negative value of \, the ws 
function embodies the following properties (see Fig. 
() a) w Q, at points of the variable-thrust sub 
irc, defined by Eq. (54); (b) wx > O, in the region A 
bounded between the m axis and the subare defined 
by Eq. (5+); and (¢) ws 0, in the region B bounded 
between the I” axis and the subare defined by Eq. (54 

On the other hand, for a given positive value of X, 
the signs which the ws function takes into the two 
regions A and B are reversed with respect to those 
indicated in Fig. 6. 

8.1.2) Analytical Proo 


The first step for carrying 


out the sufficieney proof is to analyze the class of ad 


missible displacements [Eq. (3S) and inequalities (29 
ind (39)| of the rocket in the mass-velocity plane. 
Considerations analogous to those of Section (7.1.2 
must be carried out. By drawing the four limiting 


lines 6 Q and @ B starting at | and arriving at 


2 the broken line 1362451 1s determined: it surrounds 


i closed region to which all admissible paths belong (see 


The second step is to analyze all possible combina 
tions of subares. In this connection, the analysis shows 
that the optimum path 1562 includes a subare 15, a 
subare 56, and subare 62. The equations of these 
subares are the same as in Table 1, the only difference 
being that w is replaced by wx. 

Assume now that the five quantities 77, 1), m2, V2, h 
ire prescribed and that the particular extremal path 
satistying the above conditions has been found. As 
sume also that the value of X associated with this par 
ticular extremal path is negative and that the boundary 
conditions are such that the Eulerian arc is of the tvpe 
IIT oi Table 1. 

Consider now the arbitrary, though admissible, com- 


parison path 172, and form the difference 


mY f me 


} Psd 1° + VWedm) + 


By means of Green's theorem, the above cyclic integrals 





(S4 
are trans'ormed as follows: 


Ar = {| we dVdm — l{ wad dm (S5 
SB Jd SA 


Je 
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In view of the fact that wy 0 in the region S, and ws 


> 0 in the region S,, it is concluded that 
Ar = At + XAh < 0 (S6 


The isoperimetric condition concerning the altitude rise 
requires that the contro] path 172 be chosen in such a 


way that 
Ah Ihisee hy: 0 87 
As a consequence, the inequality (S6) leads to 
At < 0 QR 


The following theorem is, therefore, proved: “All 
Eulerian ares 1562 such that the variable-thrust subare 
56 is associated with a negative value for \ minimize the 
time of flight for given altitude rise. 

By means of similar techniques it can be proved 
that: “All Eulerian ares 1562 such that the variable 
thrust subare 56 is associated with a positive value for 


\ maximize the time of flight for given altitude rise. 


CONCLUSIONS 


The indirect methods of the calculus of variations 
are employed in the study of linear problems of both 
the simplest type and the isoperimetric tvpe, subject 
to additional constraints involving inequalities. Dis 
continuous solutions are found. The Legendre-Clebsch 
condition fails to supply any information for all sub 
ares forming one extremal arc. Moreover, the Weier 
strassian function is zero along a major portion of the 
Eulerian are. 

By the use of a new and hitherto unexploited tech 
nique, based on Green's theorem the afore-mentioned 
difficulties are overcome. Rigorous sufficiency proofs 
are developed 

It is to be observed that this new technique is oi 
limited interest as far as the calculus of variations is 
concerned, because it only applies to linear problems 
Nevertheless, as far as the mechanics of flight is con 
cerned, the use of Green's theorem is rather important 
in view of the frequent occurrence of linear problems 
in the path analysis of turbojet and rocket-powered 


vehicles. 


REMARKS 


In a recent paper'* a special transformation of co 
ordinates has been indicated such that the Legendre 
condition yields a result other than zero along one of 
the subares forming the extremal arc. 

In a further paper!‘ this writer has stressed that, in 
connection with the level flight of a rocket-powered 
aircraft, Green's theorem supplies a complete proof of 
the necessary and sufficient conditions for the ex 
tremum; while the Legendre condition a/one is only 
indicative of the /oca/ behavior of the solution with 
respect to weak variations. This writer has also shown 
that, for level flight, the results of reference 13 can be 
obtained as a by-product of the approach based on 


Green's theorem 





590 JOURNAL OF THE AERO 


REFERENCES 


! Bliss, G. A., Lectures on the Calculus of Variations; The Uni- 
versity of Chicago Press, Chicago, 1946. 

2 Cicala, P., and Miele, A., Generalized Theory of the Optimum 
Thrust Programming for the Level Flight of a Rocket-Powered Air- 
craft, Jet Propulsion, Vol. 26, No. 6, pp. 443-455, June, 1956. 

3 Miele, A., An Extension of the Theory of the Optimum Burning 
Program for the Level Flight of a Rocket-Powered Aircraft, Journal 
of the Aeronautical Sciences, Vol. 24, No. 12, pp. 874-884, 
December, 1957. 

* Miele, A., and Cavoti, C. R., Generalized Variational Ap- 
proach to the Optimum Thrust Programming for the Vertical Flight 
of a Rocket, Part II, Application of Green's Theorem to the Devel- 
opment of Sufficiency Proofs for Particular Classes of Solutions, 
ZFW, Vol. 6, No. 4, pp. 102-109, April, 1958 

5 Miele, A., Soluzsioni Generali di Problemi di Ottimo in Volo 
Non-Stazionario, L’Aerotecnica, Vol. 32, No. 3, pp. 135-142, 
1952 (Translated as NACA TM 1388, 1955). 

® Miele, A., Problemi di Minimo Tempo nel Volo Non-Stazion 
artvo degli Aeroplani, Atti della Accademia delle Scienze di Torino, 
Classe di Science Fisiche, Matematiche e Naturali, Vol. 85, 
pp. 41-52, 1950-1951. 

7 Miele, A., and Cappellari, J. O., Jr., Approximate Solutions 


SPACE 


SCIENCES—SEPTEMBER, 1958 


to Optimum Flight Trajectories for a Turbo-Jet Powered Aircry 


Purdue University, School of Aeronautical Engineering, Repo 


No. A-57-7, December, 1957 

* Fraeijs de Veubeke, B., Trajectoires Optimales des Fusé, 
Revue Générale des Sciences Appliquées, Vol. 3, No. 2, pp. 33-4) 
1956. 

% Miele, A., On the Brachistocronic Thrust Program for 
Rocket-Powered Missile Travelling in an Isothermal Mediny 
Jet Propulsion (in publication ) 

” Miele, A., Optimum Climbing Technique for a Rocket-Power 
Aireraft, Jet Propulsion, Vol. 25, No. 8, pp. 385-391, August 
1955. 

'' Miele, A., The Application of Green's Theorem to Variation 
Problems of Linear Type, Douglas Aircraft Company, Inc., SM 
Report (in publication ) 

' Hibbs, A. R., Optimum Burning Program for Horizonte 
Flight, Journal of the American Rocket Society, Vol. 22, No. 4 
pp. 206-212, July-August, 1952. 

13 Ross, S., Minimality for Problems in Vertical and Horizont 
Rocket Flight, Jet Propulsion, Vol. 28, No. 1, pp. 55, 56, January 
1958. 

‘4 Miele, A., Minimality for Arbitrarily Inclined Rocket Ty 
jectories, Jet Propulsion (in publication) 


Thermal Buckling of Solid Wings. . . 


(Continued from page 580) 


distribution differs materially for some of the cases 
considered. Also, the shape of the cross section is not 
very significant, as long as the wing is thin. The 
curves of Fig. S are thus a good indication of the end 
effect in general. It should be pointed out that for 
wings of parabolic or wedge shaped section the curves 
are not applicable to aspect ratios below 2.5, since then 
the torsional analysis ceases to be valid. 

The results obtained indicate that the effect is three 
to five times as large as that mentioned by Budiansky 
and Mayers.! The error in neglecting the end effect 
can therefore be significant even for larger aspect ratios, 
and the end effect factor should therefore be included 
also in rough calculations. However, since \ can be 
taken from Fig. 8 in most cases, no additional work re- 


sults. 
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Shear Flow of a Viscoelastic Fluid Past a Flat 
Plate With Suction* 


Anadi Shanker Gupta 


Department of Applied Mathematics, Indian Institute of Technology 
Kharagpur, India 


April 15, 1958 


"| SHE PROBLEM of the laminar boundary laver on a semi-infinite 
flat plate in a uniform shear flow —i.e., in a stream of con 
Recently Sakurai 


has extended the problem to the case of an infinite flat plate with 


stant vorticity — has been investigated by Li 


uniform suction. The purpose of this note is to give an exact 
solution of the problem considered by Sakurai in the case of a 
viscoelastic fluid 

The full stress tensor p;; at any point of a viscoelastic medium 
or Maxwell body in the terminology of Reiner?) is considered to 
be decomposed into any isotropic part —pé,;; associated with iso 


tropic dilatation (actually ignored in the incompressible case 


and an additional partial stress ¢;, associated with distortion at 


a constant volume. The partial stress satisfies 


oi, + r(D/Dio:; = Qype 1 /3)M45 1) 
where 7 is the elastic constant of the medium (having the dimen 
sion of time) and the term « ] Aj,, on the right-hand side 
of Eq. (1) ts the partial rate-of-strain tensor representing the dis 


tortion only, A is the dilatation, and D/ Dt denotes total deriva 


tive with respect to time \lso ¢,, is given by 


1/2) |[(Ov,/Ox + (Ov, /OX 2 


( = 


Taking the Y-axis along the plate and the }-axis perpendicular 
to it, the equations of momentum for a steady two-dimensional 


incompressible flow of a viscoelastic fluid are 


* 1 would like to thank Dr G n the 


ration of this note 


Bandyopadhyay for his help prepa 
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p{u(du/Ox) + v(Ou/dy)] = (Op/Ox) + (00,,/Ox) + 


(O¢,,/Oy¥) (3) 


p|u(ov/Ox) + v(O0v/Oy)| = 


-(Op/Oy) + (00,,/0x) + (O0,,/Oy) (4) 


where p is the pressure and the components o-,, ¢ 
a,, satisfy (1). 


The equation of continuity is 


(Ou/Ox) + 


For the flow past an infinite flat plate, conditions will depend on y 


only. Henee Eqs. (3)-(5) become 
pu(du/dy) = —(Op/Ox) + (da,,,/dy) 6) 
pu(dv/dy) = —(Op/dy) + (da,,,/dy) (7 
(dv/dy) = 0 (8) 
respectively. Eq. (8) gives 
v = vw = constant (9) 
Hence Eq. (7) becomes 
0 = —(Op/dyv) + (do,,,/dy) (10) 
Eq. (1) now gives by virtue of Eq. (9), 
Orr + tr%(do,,/dy) = OV 11) 
or, + ru(do,,/dy) = p(du/dy) (12) 
Gy, + tw(do,,/dy) = 0 (13) 


Taking o,, 
and this case seems to be the physically interesting one 
now implies 0p/Oy = 0 while Eq. (6) shows that 0p/p-x is a func 


tion of yonly. Thus 0p/dx should be a constant (a,say). Hence 
Pp = ax + constant (14) 
Using Eq. (14) and v = % Eq. (6) becomes 
pvy(du/dy) = —a + (do,,/dy) (15) 
Elimination of o,, between Eqs. (15) and (12) gives 
(privy? — pw)(d2u/dy?) + pw(du/dy) + a = 0 (16) 
Putting r = 0, Eq. (16) reduces to Eq. (8) of Sakurai’s paper. 


The solution of Eq. (16) subject to the boundary condition 


u=Q0 aty = 0 (17) 
is found to be 
“u=l + Woy — Uo exp [voy/(v — rio?) (18) 
where (’. is an arbitrary constant and 
Wao = —(a/ pt) (19) 


If it is assumed that v < 0 so that there is a uniform suction 
at the plate, then the last term of Eq. (18) shows the exponential 
decay as y increases—viz., 


+ Woy (20) 


u—> Ua 


provided v > rz?. This indicates that the solution obtained in 
Eq. (18) corresponds to the boundary layer over a flat plate in a 
shear flow defined by 

U = Ve + Woy (21) 
the thickness of the boundary layer being of the order of (vy — 
TU") /v. Thus it is seen that the effect of elasticity in a shear flow 
is to decrease the thickness of the boundary layer. 

It is worth noting that the relation (19) among the pressure 
gradient a, the vorticity wo, and the suction velocity v is the 
same as that of Sakurai despite the effect of elasticity of the 
liquid. It may also be pointed out that the boundary-layer na- 
ture in this problem may also appear in the case of fluid injection 
v% > 0 provided vy < rz? as is evident from Eq. (18). This fact 
seems to be a departure from Sakurai’s result—viz., boundary- 
layer nature appears in the case of suction and not in the case of 


fluid injection. 


SPACE 


», and @,, of 


(Ov/Oy) = 0 (5) 


= ¢,, = 0, Eqs. (11) and (13) are identically satisfied 
Eq. (10) 
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The Effect of Extreme Cooling and Local 
Conditions on Boundary-Layer Transition | 


John R. Jack and Richard J. Wisniewski 
Lewis Flight Propulsion Laboratory, NACA, Cleveland, Of 
February 24, 1958 


[ A RECENT NOTE in the Readers’ Forum of the JOURNAL, va 
Driest and McCauley! present data showing that, under co; 
ditions of extreme surface cooling, the flow on a 10° include 


angle cone is completely laminar. This set of data appears t 


be contrary to the results presented in references 2 and 3, 


wher 
it was shown that early transition accompanied large amounts 
of surface cooling 

The purpose of the present note is to reemphasize the fact that 
mean. that 


extreme surface cooling alone does not necessarily 


early transition will be observed. On the contrary, transitior 
reversal is obtained only under certain conditions with surfac 
cooling. Reversal, as caused by single roughness elements, w 
reference 4. The 


that 


first reported by van Driest and Boison in 


» 


studies of references 2 and transition re 


influ 


3 demonstrated 
versal also occurs on very smooth bodies and is greatly 
enced by surface roughness, pressure gradients, and unit Reyn 
olds Number 

At present, all the conditions necessary for early transition ar 


not known However, a plausible explanation for the apparent 


contradiction of the reversal data of references 2 and 3 witl 
the data of reference 1 is possible if the two sets of data ar 


compared 


Before the data are compared, some of the conditions afiec 
ing transition reversal will be illustrated with data obtained from 
2°)° included angle cones (Fig. 1). The surface 
than 12 


blunted 10° ane 
microin. and_ the 
blunted tip diameter of each model was 3/16 in. Plotted in Fig 
1 is the wall to temperature ratt 


against the local unit Reynolds Number per foot 


roughness of each model was less 


average local free-stream 


The solid 
symbols designate temperature ratios at which completely lami 
nar flow was obtained, whereas the open symbols denote th 
which transition occurred for 


average temperature ratio at 


given unit Reynolds Number. Thus, for the same temperature 
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I 2 Effect of local conditions on boundary-layer transition 
on sharp cones 
tio, reversal may or may not be present, depending on the unit 
Revnolds Number and local Mach Number 


Experimental data obtained from the sharp cones are pre 
sented in Fig. 2 \lthough the trends here are not so well de 
fined as for the blunt models, a Mach Number effect is indicated 
Also included for comparison purposes is the data point from 
reference 1 Now, by considering the combined effect of unit 
Revnolds Number and Mach Number, it appears that the data 
point of reference 1 does not contradict the results of references 
2 and 3 lo be more specific, if the investigation of reference 

had been conducted at a higher unit Reynolds Number or a 
wer temperature ratio, it appears that early transition would 

ive been observed 

\lthough it is certain that the unit Reynolds Number effect 
nust be considered in conjunction with the temperature ratio 
when transition reversal is discussed, the absolute magnitude 
f this effect is not well defined. In particular, it is dangerous 
to compare results from different facilities since the unit Reyn 
Ids Number is a dimensional number. Therefore, it is not 
clear just how much higher the unit Reynolds Number or lower 
he temperature ratio of the experiments reported in reference 1 
would have to be before reversal is encountered 

In summary, extreme cooling alone does not necessarily guar 
intee transition reversal, but in conjunction with other param 
eters, such as unit Reynolds Number, Mach Number, surface 
roughness, and pressure gradients, very early transition may be 


ibserved 
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W' have recently made further observations with our smooth 
10° cone at a Reynolds Number per foot of 10 X 108 at 
local Mach Number 2.71 and still found no reversal. The wall 
to local stream temperature ratio was drawn down to 0.60. 
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Unfortunately, this Reynolds Number per foot was the maximum 
for the 12-in. supersonic wind tunnel of JPL in which we were 
testing 

This no-reversal observation was taken before ice appeared on 
the model However, as ice began to build up, transition flickered 
back and forth along the model 

We are still of the opinion that either surface irregularities or 
other extraneous disturbances are the prime cause of transition 


reversal, when it occurs, under conditions of extreme cooling 


Determination of Distribution of Twist of a 
Straight Wing to Correspond to the 
Aerodynamic-Load Distribution on a 
Sweptback Wing 


S. N. Chaudhuri and K. S. Nagaraja 

Department of Aeronautical Engineering, Indian Institute of Science 
Bangalore, India 

April 28, 1958 


| HE INTEGRAL EQUATION 


aft | 
J Vo = (1/22 Vo) § | ~,(K’) |dk K— K + a5, (kK); ] 


determines the load distribution of a sweptback wing. The 
additional downwash term o7,(«) takes into account the influence 


of the kink at the central section. The general expression 


LN = ¢ ] =a K 


for the chordwise load distribution gives the expressions for the 
effective angles of incidence in the case of finite straight and swept 


back wings, respectively, as follows 


G. = €,/4 = (c;/4xN) (1 — xh cot eM 2 
Qe = ( a ( sin ry /4arn cos ¢ * 
1 + wn(cot my — cot mH 3 


where 7, mo, and m are the values of the shape parameter », for 
straight finite wing, sweptback wing of large aspect ratio, and 
finite sweptback wing, respectively 

A sweptback wing, tapered or untape red, is sheared into a 
straight wing. Span, chord distribution, and aspect ratio re 
main unaltered 
the root and an increase near the tip, each section of the sheared 


Since sweepback causes a decrease in lift near 


wing is given a suitable twist to maintain the same spanwis« 


loading on both the sweptback and the sheared (straight) wing 
The relative angle of change in the effective incidence is given by 
6 = &./a, = [(1 — mn cot mn)/n| X 

\n cos ¢ge/{1 + mn(cot mmo — Cot mm } 1/sin mn } 


In order to apply this twist distribution 6 in the integral equation 


+1 
(2b/aye) y(n) = a — (1/2e) ff x 
] 


[d+(n’)/dn’| [dn'/(n — n’ 5 


6 


where the subscripts ¢ and 0 refer respectively to swept and un 
swept wings, and the subscript / indicates the induced incidence 


The solution 


m 


a, = —b,y, 4 > Duns 
n=1 


of Eq. (5), written in the matrix notation, gives for both swept 


and straight wings 


[Ay] yt = } aot and [A yi = ) aos 
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Since the number of 


ly} being the same in both the cases 


spanwise stations chosen and the chord distribution are the same 
in both the cases, we have 


[Ay] = |Ao] = [A] 


¥ 


ta} = [A]! fag} = [A]! feof 


Because we assume that each element of the column matrix | ag} 
is unity, we have from Eq. (6) }y} = [A]~! {6}. From y = 
cc, /2b, the distribution of local lift can readily be obtained. 

This method has been applied to a wing analyzed by Multhopp 
by his lifting-surface method,? and the results are compared 
(Fig. 1). The same order of accuracy is gained with much re- 
duced computing hours. The method has been applied for an- 
alyzing wings of various sweepback angles, aspect ratios, and 
taper ratios as indicated below: 

er.a:: 40°, 45°, 30°. 56°, 60° 
A: 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0 
A: 0, 0.25, 0.50, 0.75, 1.0 


The twist distributions calculated in these cases are plotted, 
and one carpet is given (Fig. 2) to indicate the variation of twist, 
taking the angle of attack at the central section to be 1°. For 
other angles of attack, the twist varies linearly. The carpets 
drawn will enable the determination, for practical purposes, of 
the distribution of twist for wings whose geometry falls any- 
where in the range of the A, g,.¢. and A indicated above. Since 
the values of 6 behave regularly, extrapolation is possible to a 
few values of g,.¢. and A on either side of the limits indicated 
in the table above. 
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Differential Equations for Cylindrical Shells 
With Arbitrary Temperature Distribution 


P. P. Bijlaard* 

Consultant, Bell Aircraft Corporation, Buffalo N.Y., and Professe 
Cornell University, Ithaca, N.Y 

April 21, 1958 


— TEMPERATURE is assumed to vary linearly in the radj 
direction Let 7 be the temperature rise in the middle play 
and let the outside temperature be 7) higher than the inside ten 
perature, both 7 and 7) varving arbitrarily in the longitudin 
Y-) and circumferential (¢-) directions 

With a coefficient of thermal expansion a, 7 causes strains a] 
in all directions. Hence the stress causing strains in the middk 
plane are «, — a7 and e — a7, where e& and e: are the tot 
strains in the middle plane in the Y- and ¢-directions. Th 
temperature difference 7; causes curvatures a7)/h of the she 
wall, with thickness 4, so that the stress causing curvatures ar 
Xx — (a@7\/h) and xg — (a@7\/h), where x, and xg are the tota 
curvatures in the Y- and g-directions. Hence in the formulas 
for the membrane forces, .V, and Ny and the bending moments 


¢¥ 





M, and Mg, 4, €, xz, and x, have to be replaced by « — a7 
e — al, xr — (aT\/h), and xg — (a@7\/h), so that 
N, — Eh (1 — v?)} le ~~ tse = (l + pv al 
No = |Eh (1 — v®)| [es + vey — (1 + vial ; 
M, = —D[xz + vxg — (1 + v)aTy/h] 
My, = —Dixe + vxr — (1 + vali /h 


where D = Eh'/{12(1 — v?)| is the flexural rigidity and » is 
Np, Wy, Miyy, and My, do not change. Ex 
pressing €,, €, x2, and xy in terms of the displacements uw, 7, and u 


Poisson's ratio exe 


in the Y-, g- and inward radial (Z-) directions, Eqs. (1) ean b 
written as 


N, = F, — EhaT/(1 — » | 
= F, — EhaT/(1 — » 

M, = F; + (1 + »)Dal; | 

Mo. = Fy + (1 + v)DaT\/h 


where F, to F; are the same functions of u, 7, and w that appear 
in the formulas applying for zero temperature rise. From refer- 
ence 1, the three equations of equilibrium are 


aN.’ + Ngx° + aX =0 } 
an,” + PN - M,° —a Me! + a*> = 0) 3 
My? + aMygy’? + aMgx"? + a2M,” + aN, + a2Z = 0\ 


where Y, ®, and Z are the loadings per unit surface in the X-, 
g-, and Z-directions, a is the shell radius, and the symbols 

and (°) indicate partial differentiation with respect to x and ¢ 
If X, ®, and Z are assumed to be zero and Eqs 
(2) and the unchanged expressions for Nyy, Nox, Wex, and Mos 


respectively. 


are inserted in Eqs. (3), one obtains 


F; — |Eha/(1 — v)laT’ = 0 | 

F, — |Eha/(1 — v)Jal[T° + (ka/h)T,°] = 0 4 

F; + [C1 + v)Da/h] (a?T,” + 7™°°) - | 
[Eha/(1 — v)laT = 0 


where F;, F;, and F; are the same functions of 4, v7, and w that ap 
pear in Flugge’s three simultaneous partial differential equations, 
applying for zero temperature rise and k = h?/(12 a”). On the 
other hand, for zero temperature rise and Y, #, and Z unequal to 
zero, Eqs. (3) become 


F, + aX = 0 / 
Fe + a? = (> 5 
F:+ az = 0\ 


where F;, Fe, and F; are the same functions as those in Eqs. (4 
From the definitions given in the foregoing 


* The writer wishes to thank the Bell Aircraft Corporation for its permis 
sion to publish this note, and Mr. Arthur Schnitt, Asst. Manager, Develop 
ment Analysis Department, Special Weapons Division, for having asked 
him to consider this problem 
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Eh (1 — v) = (1 + vw) D/(atk 
that from comparison of Eqs. (4) and (5) the effect of the tem 
nerature rises 7’ and 7) is equivalent to that of surface loads 
\ - 1 + v)D/(a*k)|aT”’ | 
>, = 1 + v)D/(a*k)|al{7 Tt ka /h)T, 6 


- 1+ Pp Dita k)|al|7T — (ka h a?7,” + 7; \ 


Using these fictitious loads in the Donnell-type equations de 
rived in reference 2, considering Rk as small with respect to one 
ut not neglecting terms where & is multiplied with derivatives of 


her order than that of terms without &, one obtains the fol 


ng 


wing equations 


—" 2varu + baiw’’’ + 204 — viarw’ + 
2 12(1 — v?)/h®la'w’’’’ + (4 — Br? )atw'’”’ 4 
2(2— yp aru ‘ T « T 1 + viaa k x 
PviT — (kab /h)V'T, + T°? — va?T" + 
hor. 7" =0 (7a 
v ~ + kala -— ow + 
1+» aaja?l 7 7 - 
1+ v)/(1 — v)] (ka/h)7 { (7b 
v 2 — viarx + WwW — 2ka*(a*w'’’’? 4 
ww’ T 1+ p aa\a?] 9 + 7 T 
2 l1— » ka*/h)T,"’ +- (ka/h)7 i iC 
where v2" = (0?/Ox? + 0?/a20¢?)’ 


miditions, 


Depending on the type of loading and boundary c 
some of the terms in Eqs. (7) may become unimportant and can 
then be dropped 

If the shell is simply supported at the ends and 7 is the tem 
perature rise above that of the heads, which is assumed to be 
constant, Eqs. (7) can be solved by expressing T, 7), w, u, and ¢ 
is double Fourier series. This leads to similar results as derived 
in reference 3 where in the present case, of course, the extra 
terms shown in Eqs. (2) have to be added. Also for clamped 
shells, a solution was worked out using the same normal functions 
(, that were used in reference 4 for determining the buckling 


load of clamped shells 
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Flutter of a Thin Membrane in Hypersonic 
Flow 


Franklin D. Hains 

Aerodynamicist. Research Division, Bell Aircraft Corporation, 
Buffalo, N.Y 

May 7, 1958 


pepe A two-dimensional membrane with one side exposed 
to the air stream while the other side is kept at free stream 
pressure (see Fig. 1 This is the limiting case of a thin panel 
with zero bending stiffness or, under certain conditions, might 
be the membrane-like hypersonic sail proposed by Daskin and 
Feldman The differential equation governing the small vibra 


tions of such a membrane is 


MS — {srr = PD. —P (1 
where f is the tension, m is the mass per unit area, / is the time, 
ind p p., is the net aerodynamic pressure in the 2 direction 


a>, and if w/ax < 1 so 


If the flow is hypersonic, i-e., 
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Fic. 2. Characteristics network 


that the flow field is free of shocks, the pressure given by piston 


theory? is 


Pu — Po = Pata? (w/a + l(y + 1)/4](w/a,)? 4 
(y + 1)/12|(w/a ; (2 
where w = 2, + uyzr. Eq. (1) is linear if the first-order term in 


Eq. (2) is retained and the others neglected. An exact solution 
was given by Goland and Luke,’ who pointed out that lengthy 
numerical calculations were still necessary in order to study a 
particular problem. If the second- or third-order terms are re 
tained, Eq. (1) is nonlinear and exact solutions are not readily 
obtained 

The purpose of this note is to present a straightforward numeri 
cal method for obtaining solutions of Eq. (1) when one or mor 
terms in Eq. (2) are retained. It is based on the method of 
characteristics,‘ and, in the opinion of the writer, gives more 
physical insight into the nature of the oscillations 

By introducing the variables 


Zesz/L, X=x/L, T= a,t/L 
V = m/p UB F = f/p_a ef... P= p pa 2) 


Eq. (1) is transformed into the nondimensional form 


MZrr — FZxx =P, —P (4 
while Eq. (2) becomes 
P, — Po =r+ (y + D/A? + My + ODD/12) (5 
where 
d g+uh/a,. ¢£ Zr, Ah Zx 6 


Eq. (4) is hyperbolic and the characteristic lines in the Y-7 plane 
are given by 

dX/dT = +V F/M (7 
where the positive sign designates the / family. The corre 


sponding compatibility relations are 
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dh + |(P P,)/F\dX = V M/F dg = 0 (8) 

The characteristic lines, whose slopes are known, are straight 
lines because F and \/ are constant. The characteristic network 
in the Y-7 plane is shown in Fig. 2. Because the membrane is 
fixed along Y = OQ and Y = 1, the characteristic lines are reflected 
along these lines. As was expected from physical reasoning, the 
Dirichlet boundary condition Z(0, 7) = Z(1, T) = 0, and the 
Cauchy conditions, Z and Z; given for 7 = 0, are sufficient for 


this U-shaped boundary. 

The compatibility relations hold along the characteristic lines 
and are expressed in finite-difference form. Using first differ- 
ences and referring to the numbered grid in Fig. 2, 


/ 4 
[M/F go 
v p= (/2KHF (Pu — Pa — (Pu — Pa dsl/Fl X 
((X; — X,)/2] + WV M/Fles = 2) = hy +h) (9) 
By proceeding stepwise from 7 = 0, the values of g and h are 


found over the X-7 plane. The integral surface Z(Y, 7) is 
obtained by a simple quadrature. 
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Effect of Axial Loads on Lateral Vibrations of a 
Slender Member With Any Degree of 
End Restraint* 


G. M. Smith and L. E. Young 


Associate Professors, Engineering Mechanics, The University 
of Nebraska, Lincoln, Neb. 


May 7, 1958 
SYMBOLS 

£ = length of member 
kL = distance between inflection points 
A = time function 
V total potential energy 
7 = total kinetic energy 
E = modulus of elasticity 
I = moment of inertia about axis perpendicular to plane of vibration 
P = axial load 
P. = critical buckling load = r*EI/k?L? 
Wo = moment at supports of member 
Y = mass per unit length 

= velocity 
yo = displacement of end supports from reference axis 
p = circular frequency 
e = WEI/y 
f frequency of member 


mM” of the theory of vibrations of axially loaded members has 
been confined to the extreme end conditions—i.e., pinned 
and fixed. In reality very few axially loaded members can be 
classified as either pinned or fixed. This paper develops a rela- 
tionship between lateral vibrations as affected by axial load and 
any degree of symmetrical end restraint. 

The members considered are long, slender, axially loaded mem- 
bers with constant stiffness and symmetrical conditions of end 
restraint. Only the fundamental mode of lateral vibration is 
considered and the effect of damping is neglected. A simple 


* This investigation was sponsored by research grants from the George 


Abel Memorial Fund. 
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Fic. 1. Configuration of member with any degree of sy mmetri 
end restraint 


shape function is obtained by employing a reference axis whi 
moves with respect to the end supports as shown in Fig. | 

The positicn of the reference axis is determined by the degr 
of end restraint index, k, where k may have any value betwe 
0.5 (fixed ends) and 1.0 (pinned ends 

The solution is based on energy considerations and Lagrang: 
equation 

Fig. 1 shows the shape configuration for an axially loade 


member with end moments whose magnitudes depend upon th 


degree of restraint. The shape function is taken as 
y = A[l — cos (xx/kL)] 
In Eq. (1), -1 is a function of time, and the shape function satis 
fies the following boundary conditions with respect to the refer 
ence axes: 
y=0 atx =0 andx = 2kL 
dy/dx =0 atx =0 andx = 2kL 

d*y/dx? = M/EI=0 atx = kL/2 and x = 3kL/2 
The potential energy V of the system between supports is 
follows: 


(L/2 


PkL+(L/2 
V = (EI /2) f, (d?y/dx?)*dx + Mo(dy/dx)\x =Ki- (1/7 
L- 


From Eq. (1), 
dy/dx = A(m/RL) sin (wx/kL) 3 
d?y/dx? = M/EI = (Ax?/k?L?) cos (xx/kL } 
M)/EI = d*y/dx? xy =k1n_(1/2 5 
Substituting Eqs. (3), (4), and (5) into Eq. (2) and integrating 
the potential energy becomes 
V = A? 4} (EI 13/2 R3L3) [(4/2k) — (1/2) sin (2/k)] — 
(Px/2kL) [(w/2k) — (1/2) sin (r/k)]{ (6 
The kinetic energy of an infinitesimal length dx is (1/2) 70d 
where v is the velocity of the elemental length with respect t 
an axis through the ends of the member. 
v = dy/dt — dy)/dt\x =k 
Then the total kinetic energy is 
: kL+(L/2 
T = (3 ad? f, [cos? (4/2k) + 2 cos (4/2k) X 
kL—(L/2) 


cos (rx/kL) + cos? rx /kL\d 


where denotes dy/dt. Integrating 


T = A2(yL/2) [cos? (r/2k) — (2k/x) sin (4/k) + 
(1/2) + (k/2m) sin (/k 5 


Using the potential energy from Eq. (6) and the kinetic energ) 
from Eq. (8) to satisfy Lagrange’s equation, 


(d/dt)(OT/0A) — (OT/0A) + (OV/OA) = Q 9 
The following ordinary linear differential equation results: 
EIr | (w/2k) — (1/2) sin (2/k) ( 


yk8L4 V1 + (1/2) cos (x/k) — (BR/2m) sin (4/k) § 
[1 — (P/P,)JA_ (1 


x 
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Fic. 2. The relation between load and natural frequency for 
different degrees of symmetrical end restraint 
Satis ; 
refer The solution of Eq. (10) is 
A = Csin pt + D cos pt 
where 
: Elx* 5 (w/2k) — (1/2) sin (2/k ly 
p= L 4 
2 >k8L41 + (1/2) cos (r/k) — (3k/2m) sin (2/k { 
is 1 — P/P.]) 
Then the frequency becomes 
2 f= P/2xe = g/2kL? X 
j (w/2k) [(a/k) — sin (2/k)] (2 
2 a. 
) 1 + (1/2) cos(4/k) — (3R/2m) sin (x/k) § 
(1 — (P/P.)]”? (11) 
Eq. (11) reduces to the solution for pinned ends when & equals 
4 one and for fixed ends when & equals one half. Fig. 2 is a graphi- 
5 cal representation of Eq. (11) and shows the effect of axial load 
: and end restraint on the frequency for any member of given 
ing . . . 
length, modulus of elasticity, moment of inertia, and mass per 
unit length 
The authors believe that with additional analysis for nonsym- 
6 metrical vibrations and with careful technique it would be pos- 
rs sible to determine the load on a column in an existing structure 
2d , - : a 
from the vibrational characteristics of the column 
t 
- + 
Reynolds Number Effect on Small-Engine 
Compressor Performance 
d. 
Frank E. Lenherr 
Small Aircraft Engine Department, General Electric Company, 
Lynn, Mass. 
x May 12, 1958 
2 [' IS WELL-KNOWN that air-breathing propulsion systems oper- 
ate differently at altitude conditions than at sea-level static 
conditions. This effect has been referred to as altitude effects! or 


Reynolds Number effect (references 2-4, for example) depending 
upon the particular phenomena being investigated. 

The axial flow compressor has been the turbjojet component 
that has received the most critical scrutiny since it is the most 
complicated and difficult component to comprehend. However, 
no satisfactory theoretical method currently exists for predicting 
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the altitude effects on compressor performance prior to test 
evaluation. All valid predictions rely upon previously obtained 
test results from similar designs 

The advent of the smal!, high-performance, turbojet engine 
posed a serious question regarding the effect of reduced densities 
In order to actually determine the effects, a series of tests were 
conducted on multistage engine compressors of small engines 
The variation in chord length versus axial location 1s shown in 
Fig. 1 for three engine compressors tested by this department 


The measured variation in overall compressor efficiency in 


terms of Reynolds Number index, (Re), je /( Re dee el (Se€ 
Fig. 2) showed only a moderate drop off in overall efficiency 
ratio, ()attitude/(1's compared to previously published 
data.? 


These tests have conclusively demonstrated that small engines 
can be designed to vield good altitude performance without re 
sorting to the use of large chord lengths 

Significantly, these data tend to confirm the previously pro 
posed hypothesis that the type of stage loading distribution, 
front-stage design, and turbulence level® in the compressor are 
undoubtedly the controlling factors. Further investigations are 
required to evolve a satisfactory theory for altitude performance 


prediction 
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A Numerical Method for Solving Boundary- 
Layer Equations 
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Graduate Schoo! of Aeronautical Engineering 
Cornell University, Ithaca, N.Y 


May 19, 1958 


M”: BOUNDARY-LAYER PROBLEMS lead to differential equa- 
tions whose solutions can only be obtained by numerical 
means. Through the years, a number of methods have been 
developed; each, no doubt, has its value, and a comprehensive 
study of their relative merits is a task that remains to be under- 
taken by numerical analysts. In the present note, we wish to 
outline a method which has been used with considerable success 
in obtaining numerical solutions to equations commonly encoun 
tered in boundary-layer problems, and which is particularly 
adaptable to the high-speed computers currently in use 

The differential equations resulting from most boundary-layer 
investigations are ordinary linear equations of the third order 
with known variable coefficients, the general form of which may 
be written as follows: 


+ agf’’ + bg’f’ +ce"f =H 


where f represents the function to be determined, g and H are 
known functions given in tabulated values, while a, 4, and ¢ are 
constants 

The numerical solutions of this type of equation can be ob- 
tained through the collaboration of two well-known numerical 
methods: (1) Milne’s 4-Point Method, and (2) the Trapezoidal 
Formula.' * 

Although the Milne 4-Point Method generally permits greater 
accuracy, the actual application to the boundary-layer equations 
through the entire range of the independent variable ¢ is ham- 
pered by the emergence of divergent oscillations at large values of 
¢. This phenomenon is usually preceded by a noticeable increase 
in truncation errors. Investigation of the asymptotic behavior 
of the functions at large values of ¢ reveals that in general the 
asymptotic solutions have two branches, one of which behaves 
like a polynomial in ¢. This branch is normally discarded be- 
cause of finite boundary conditions at infinity. But its presence 
cannot be precluded in the process of computation before adjust- 
ment is made to fulfill the boundary conditions at infinity. This 
seems to contribute to the inherent instability of the functions 
at large values which has rendered the method unsuitable in 
that range 

Fortunately, for the equations concerned, the values of the 
functions and their derivatives usually become comparatively 
quiescent before the above-mentioned phenomenon sets in. <A 
switch is therefore devised, at a suitable point, to the Trape- 
zoidal Scheme which permits less truncation errors without im- 
pairing overall accuracy. The criterion for the crossover is that 
the overall truncation error through the entire range of computa- 
tion bea minimum. Sucha point obviously cannot be set a priori 
but can readily be located after a few trials. 

The start of the 4-point method requires four consecutive 
known values of the function f and its derivatives. These are 
usually given partially by the boundary conditions at the starting 
point and through the differential equations themselves, and the 
rest mav be calculated by Taylor’s expansion or other starting 


devices. Denoting these values by fi, fo™, f;°, and fy, then 


f;™ can first be predicted by the formula: 


$h/3) (2f@tV — fynty + 
2fo™ +1 ) + (28 90) hof(n +5)(f) 


fe = fi + 


(1) 


SPACE 
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then corrected by: 
hf” = fe + (h/3) (fs"79 + 4f,T) + fat) — 
(1/95) Hef" 75 (2) (9 


where h is the step of the variable ¢, and ¢ lies in the intery, 
spanned by each formula. The values thus obtained are fed int 
"is ob 
tained. New /;’’, f;’, f; are in turn calculated from Eq. (2) an 
the process of iteration is repeated until the required accuracy j 


the differential equation, whereby a new and better / 


obtained. ‘ 
In using the Trapezoidal Scheme, the following formulas ay 


used as predictor and corrector, respectively : 


A = fs + QA + (3/3) fe 79 (0) 


fx = fy + (h/2) (fp t) 4+ AY) — (h3/12)f +92) (4 
The process of iteration can be carried on much the same 
before 

Of the three boundary conditions generally required to soly 
the boundary-layer equations, two are normally given at ¢ =| 


and the third at ¢ —~ ©, such as: 


f(O) = 0, (0) = 0, f’( ©) = constant 


The starting values of the unknown function therefore are vet 
not entirely determinate. The coefficients of the differential 
equation are such that f’’(O) is the only value that needs to be 


assumed to start the computation. The boundary condition at 


¢— © is met by assuming a suitable value of G’’(0) to find a com 
plementary function G and another suitable value of J'’(0) to 
find the corresponding particular integral 7. Any solution of the 
equation may be written as (AG + /) where d is an arbitrary con 
stant to be determined by the condition at ¢—~ ©. The correct 
f’’(O) is therefore }AG’’(O) + 7'’(0) |. 

Since this computing method is essentially one of iteration, any 
value not given by the above formulas can be started with some 
crude trial values like fj’ = fo’ + hfi’’. 
the onset of the program, auxiliary formulas may be used to start 


To ensure accuracy at 


the computation first at very small intervals, such as: 


predictor: fg” = fa™ + [(B — a)/2] (fat? + fet) (5 


corrector: fg" = fal” + [(8 — @)/2] (fa"*” + 


The intervals (a,8) can then be successively doubled until the in- 
terval for the tabulated coefficients is reached. For example, when 
the table of the coefficients is given in (0.1 intervals, the starting 
computation may proceed with the following intervals: 0, 0.0125, 
0.025, 0.05, 0.1, 0.2, 0.3. The values of f™ at 0, 0.1, 0.2, 0.3 
thus obtained can then be used to proceed with the methods out 
lined above. 

The truncation errors of the Milne and Trapezoidal formulas 
are indicated by the last terms of Eqs. (1), (2), (3), and (4) 
The truncation error for each step of computation can be calcu- 
lated by the differences of the predicted and corrected values 
The sum of the absolute values of these errors is therefore a safe 
estimate of the upper bound of the truncation error of the entire 
run. 

One advantage of the method is that the formulas used are 
simple and direct, with minimum of arbitrariness. The trunca 
tion errors which are readily traceable at every step make for 
convenient built-in automatic control when programed on high- 
speed computing machines. This method has been suggested by 
Prof. J. B. Rosser, Department of Mathematics, Cornell Uni- 
versity, and programed for machine computing at the Cornell 
Computing Center on the IBM Card Programmed Calculator 
and the 650 Magnetic-Drum Computor. It has been extensivel) 
tested by using known functions, and no particular difficulties 
are experienced in obtaining functions with comparable accuracy) 
to that of the coefficients in the differential equations. The 
method has since been applied to a number of boundary-layer 
investigations with satisfactory results. 
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me a 
yee SMAI L-DISTURBANCE THEORY Of boundary-layer stability ‘ 
solve has attracted serious attention as a possible explanation of Hs 
es. the ‘triggering’? mechanism for boundary-layer transition from i Ea 
a laminar to a turbulent state. At low speeds the small-dis- 
turbance theory! was conclusively supported by the classic ex- Bs, 
periments of Schubauer and Skramstad? and of Liepmann.* f 
€ yet Later, Lees and Lin‘ and Lees® extended the theory to include 
ential the effects of compressibility and heat transfer, but it was only 
to be recently that the existence of ‘‘stable’’ and ‘‘unstable’”’ ranges of 
on at disturbance wavelengths were discovered by Laufer and Vre- 
com balovich® in supersonic flow. 
)) to The present experiment’ was undertaken in order to investi- 
f the gate the stability of the laminar boundary layer in hypersonic 
con flow. More specifically, the purpose was to see if there exist Fic. 2. Variation of the amplitude of a f = 5.5 ke. per sec 
rrect small laminar fluctuations which amplify or damp as they pro- fluctuation with distance x from the plate leading edge 
gress downstream in the boundary layer, and, if so, whether 
any there exist distinct ranges of the fluctuation wavelengths and the - F ; ; 
; . ‘ : The streamwise growth or damping of the fluctuation ampli 
ome mean-flow parameters for which the fluctuations are unstable in : ‘ : 
ie tudes in the boundary layer was observed with the aid of a 
y at a manner similar to the low-speed boundary-layer case. 2 : ; , ‘ 
a : , : ‘ ee ae 7 0.0001-in.-diameter platinum-rhodium (90 per cent Pt, 10 per 
tart The experiment was carried out in the GALCIT 5- by 5-in. ‘ pits , — 
hot-wire. The wire output variation (root mean 


in- 
hen 
ing 


hypersonic wind tunnel (Leg 1), at a Mach Number 5.8. The 
model employed was a sharp-lipped, smooth steel flat plate with 
a span of 5 in. and a chord of 26 in.; the plate was thermally in- 
sulated and was aligned with the flow at zero incidence. Arti- 
ficial disturbances in the form of feeble air pulses were injected 
into the plate boundary layer through seven 0.060- by 0.028-in. 
slits drilled on the plate surface along a straight line parallel to the 
plate leading edge and about 1.6 in. behind it. The periodicity 
of these fluctuations was achieved with the ‘‘siren”’ 
shown in Fig. 1, and their frequency could be set at will between 


mechanism 


l and 40 ke. per sec. 


* The work discussed in this paper was carried out under the sponsorship 
and with the financial support of the Office, Chief of Ordnance, and the 
Office of Ordnance Research, U.S. Army (Contract No. DA-04-495-Ord-19) 
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mechanism 


cent Rh) 
square signal) was recorded at each of a series of selected narrow 
frequency bands as the wire was moved downstream (along x) 
at constant relative distance y/6 above the plate surface, where 
6 is the local thickness of the boundary layer. Minima and 
maxima in the wire output thus correspond to the ‘‘lower’’ and 
‘upper’ neutral branches, respectively, in the boundary-layer 
stability diagram. 

Fig. 2 shows a typical fluctuation amplitude 
viewed on the screen of a display-type harmonic wave-analyzer; 
a minimum in the amplitude of a 5.5 ke. per sec. fluctuation is 
seen to occur at about 4.3 in. downstream from the plate leading 


Variation as 


edge. <A series of such minima and a series of maxima were ob 
tained for various combinations of the ‘‘siren’’ frequency and the 


wind-tunnel flow parameters. These neutral-stability points are 
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Fic. 3. The neutral stability branches for J = 5.8. Open 


circles denote points on the lower neutral branch; solid circles 


denote points on the upper neutral branch 
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plotted on Fig. 3, properly nondimensionalized in the familiar 
Bv,,/l,,* — Rg scheme. 

It is apparent that the shape of the instability region at J = 
5.8 is generally preserved in comparison with the stability dia- 
gram at lower Mach Numbers, although the area it encloses 
seems to be much smaller. Both neutral branches are displaced 
to considerably lower frequencies (larger wavelengths) than the 
Schubauer-Skramstad neutral branches, and to higher Reynolds 
Numbers than the corresponding Laufer-Vrebalovich data at 1/ 
= 2.16. It is also worth noting that the amplification rates ob- 
served were generally small; for example, the total amplification 
across the instability region is of the order of 1.4 at a nondimen- 


sional frequency of 5 & 107% 
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pep DERIVED necessary conditions for minimum-drag wings 
of zero thickness when the wing plan form, free-stream Mach 


Number, and lift are prescribed. He showed that on the wing 
plan form the downwash in the combined flow field must be a 
given constant dependent on the lift. Sedney? used the calculus 
of variations to arrive at the same result. He restricted the class 
of possible variations, however, to those which did not give edge 
suction forces, assuming that the Ursell-Ward reciprocal rela- 
tions were valid only under these conditions. The latter re- 
striction was subsequently shown by Jones’ to be unnecessary. 

In the present paper we shall show that the optimum con- 
figuration is not necessarily unique when one has subsonic lead- 
ing edges with the possibility of having optimum warps with 
and without edge suction forces; only in the latter case will 
Jones’ criterion be shown to be valid for supersonic flow 

Let us first denote the local lift of the wing by P and the angle- 
of-attack distribution by a. The domain of the wing plan form 
is C + S where C is the boundary of the wing. If do and d? 
denote the differential area element of S and the differential 
length of C, then the drag D of the wing may be expressed as 


D= 3 Pads — gf. A(1)G*(1, a)dl (1) 


Here in the contour integral for the suction force at the wing 
leading and trailing edges, the function G at a point / on C is the 
limiting value of essentially the pressure P multiplied by the 
square root of the distance normal to local wing edge as a point 
/ is approached from the wing interior; it depends on the angle of 
attack of the wing in the region of influence of point /o0n C. (See 
reference 4 for more details.) Additionally, the function A(/) 
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is a given function depending on the “local sweepback”’ of C ang 
the free-stream Mach Number. 
The total lift L may be written as 


L= f. Pda (2) 


Our problem is now to minimize D keeping L constant; that is, 
we form F = D + AL where d is the constant Lagrange multp 
plier, take the first variation of F with respect to a, and set the 
resulting expression equal to zero. We then obtain 


6F = { (Péa + abP + bP )\da — 
as 


of. A(IG(, a )bGdl 
( 


where 67? and 6G are given functions of da 
For the further considerations we require the reciprocal rela. 


tion of Ursell and Ward® which states that 


J Piada (4) 


where @ is the angle of attack required in the reverse flow to ob 
tain the same pressure distribution ? on the plan form. Ags 
noted before, Eq. (4) has been shown by Jones* to be valid eveg 
with edge suction forces 

Eq. (4) is now to be inserted into Eq. (3), with the result 


f (a + @+A)6Pdo — 2 £. A(1)G(, a)é6Gdl (5) 
S ( 


There are now two ways in which Eq. (5) may be fulfilled, 
First, it is possible to have each of the integrals vanish separately, 
Since 6? and 6G are arbitrary, this implies that 

a+é&= —dA = const. on § 


G= 0 on C 


The former equation expresses the constancy of the combined 
downwash (the original Jones’ criterion ), while the latter requires 
additionally no edge suction forces; that is, no leading-edge 
suction force and the fulfillment of the Kutta condition at the 
trailing edges. The second, less obvious way in which Eq. (5) 
may be satisfied is that the two integrals while each being non- 
zero just negate each other inthe sum. This is possible since 6P 
and 6G, although arbitrary, may not be independent, both being 
functions of a. In this case, Jones’ criterion will be invalid, and 
edge suction forces will be present. 

An example fulfilling the first criterion is the case of an elliptic 
wing at supersonic speeds. Here use of Jones’ criterion yields 
a constant loading on the wing; thus no edge forces will be present 
in this case 

The flat delta wing at Mach Number 1 is an example illustrat- 
ing the second possibility. (The Mach Number 1 case is con- 
sidered here as a limiting case of supersonic linearized theory.) 
It can be shown that the angle of attack in the reversed flow can- 
not be constant and have at the same time the same pressure dis- 
tribution as the forward flow. Here, clearly, would be an opti- 
mum case with leading-edge suction force not fulfilling Jones’ 


criterion 
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